{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import scipy.io as scio\n",
    "import numpy as np\n",
    "import random\n",
    "import scipy\n",
    "import scipy.sparse as sparse\n",
    "import scipy.sparse.linalg as splin\n",
    "\n",
    "import h5py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 254,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_f(num, l0, seed=None):\n",
    "    '''l0 组件尺寸，num 组件个数'''\n",
    "    num_row = int(l/l0)\n",
    "    step = int(l0/(l/N))\n",
    "#     print(f'num_row{num_row}, num{num}, step{step}')\n",
    "    f = np.zeros((N, N))\n",
    "    if seed is not None:\n",
    "        np.random.seed(seed)\n",
    "    place = np.random.choice(num_row**2, num, replace=False)\n",
    "#     print(place)\n",
    "    for i in place:\n",
    "        row = int(i // num_row)\n",
    "        col = i %  num_row\n",
    "#         print(f'row{row}, col{col}')\n",
    "        f[row*step: (row+1)*step, col*step: (col+1)*step] = -1\n",
    "    f[[0,-1], :] = 0\n",
    "    f[:, [0,-1]] = 0\n",
    "    return f # N*N\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 64\n",
    "l = 0.1\n",
    "h = l/(N-1)\n",
    "u0 = 298/1e4 \n",
    "\n",
    "B = sparse.diags(-4*np.ones(N-2)) + sparse.diags(np.ones(N-3), 1) + sparse.diags(np.ones(N-3), -1)\n",
    "I = sparse.eye(N-2)\n",
    "II = sparse.diags(np.ones(N-3), 1) + sparse.diags(np.ones(N-3), -1)\n",
    "A = sparse.kron(I,B) + sparse.kron(II, I)\n",
    "\n",
    "# B = np.diag(-4*np.ones(N-2)) + np.diag(np.ones(N-3), 1) + np.diag(np.ones(N-3), -1)\n",
    "# I = np.eye(N-2)\n",
    "# II = np.diag(np.ones(N-3), 1) + np.diag(np.ones(N-3), -1)\n",
    "# A = np.kron(I,B) + np.kron(II, I)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<3844x3844 sparse matrix of type '<class 'numpy.float64'>'\n",
       "\twith 18972 stored elements in Compressed Sparse Row format>"
      ]
     },
     "execution_count": 135,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 256,
   "metadata": {},
   "outputs": [],
   "source": [
    "F = get_f(20, l/10, seed=10)\n",
    "f = F.copy()\n",
    "f = f[1:-1, 1:-1]*h**2\n",
    "\n",
    "\n",
    "f[0, :] -= u0\n",
    "f[-1, :] -= u0\n",
    "f[:, 0] -= u0\n",
    "f[:, -1] -= u0\n",
    "\n",
    "u = splin.spsolve(A, f.reshape(-1,1)).reshape(-1,1)\n",
    "# u = np.linalg.solve(A, f.reshape(-1,1))\n",
    "# u = np.linalg.solve(A.astype(np.float32), f.reshape(-1,1).astype(np.float32))\n",
    "U = u0 * np.ones((N, N))\n",
    "U[1:-1, 1:-1] = u.reshape(N-2, N-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 224,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.883755872115e-16"
      ]
     },
     "execution_count": 224,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(A@u - f.reshape(-1,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 221,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "200"
      ]
     },
     "execution_count": 221,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 225,
   "metadata": {},
   "outputs": [],
   "source": [
    "UU = np.zeros((1,3,N,N))\n",
    "UU[0,0,:,:] = U\n",
    "with h5py.File('f.hdf5', 'w') as file:\n",
    "    file['input'] = F.reshape(1,1,N, N)\n",
    "    file['output'] = UU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 239,
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File('f.hdf5', 'r') as file:\n",
    "    F = file['input'][:1]\n",
    "    U = file['output'][:1]\n",
    "F = F[0, 0, :, :]\n",
    "U = U[0, 0, :, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 257,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x1b41043c9c8>"
      ]
     },
     "execution_count": 257,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUcAAAD8CAYAAADkM2ZpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2df7RdVXXvP18TkrQoP0yAQQkKDGghERXhgQ7bSonQ4CuCBoVoq3RkPHxgHNbf5CnWpvQZRYsZQ7SNQPlhX4HSUoOhRX5Wi2gJDaCBEgLlyRXaEIi0vpZEYL4/1jq4su7Z666zs++954b5GeOMe/Zea6+9zo87z5pzzR8yMxzHcZzteclkT8BxHGcYceHoOI7TBxeOjuM4fXDh6DiO0wcXjo7jOH1w4eg4jtMHF46O44wrkhZKekDSRknn9GmfKemq2P59SQfE88dLukvSD+Lf45JrTpN0r6T1kj6fnH+lpJtj222S5iZtr5D0LUn3S7qvd58mXDg6jjNuSJoGXAicCMwDFkual3VbAmwxs4OBC4DPxfObgZPM7HDgvcAVcczZwPnAAjObD+wjaUG85gvA5Wb2amA58NnkPpcD55vZYcDRwKbS3HdIOI71i+A4zoueo4GNZvawmW0DrgROzvqcDFwWn18DLJAkM1tnZo/F8+uBWZJmAgcBG8zsidh2E7AoPp8H3Byf39q7VxTI083sRgAz+6mZ/Wdp4tMHf62B5BfheGAEuFPSajO7r3CNh+M4zjhjZtqR6xcuXGibN2+u6nvXXXfdYGYLC132Ax5NjkeAY5r6mNmzkp4GZhNWjj0WAevMbKukjcChUS0eAU4BZsR+98S+K4G3AS+LK81fBn4i6a+BAwkC9Rwze65p4q2FI8kvAoCk3i9Co3B0HGf42bx5M2vXrq3qK+lQSWnnVWa2Ku3S57J8kVTsI2k+QdU+AcDMtkg6C7gKeB74LmE1CfBR4MuSzgC+DfwYeJYg634NOAL4Ubz2DODipte2I8Kx5hcBSWcCZ+7AfRzHmWAGyLmw2cyOKrSPAPsnx3OBxxr6jEiaDuwOPAUQN1SuBd5jZg8l87sOuC72ORN4Lp5/DHh7PP9SYJGZPS1phLDy7C3m/gZ4PeMkHGt+EYi/IqvihNJfg+0HS45f8pLtTaHpcfp8+vTpVf3GGr+pX2mONed3pG/65cy/qM8991xjW9OXuvRlb9vWpt94jNHmuq4SrnTxfjd91nm/559/vmqMLkjvtYPcCRwi6UDCKu504F1Zn9WEDZc7gFOBW8zMJO0BrAGWmdnt6QWS9jazTZL2BM4G3hnPzwGeMrPngWXAJck89pS0V7RVHgcUl8c7siFT84vgOM4Uw8yqHxVjPQssBW4A7geuNrP1kpZLemvsdjEwO9oSPwz0NneXAgcD50q6Oz72jm0rJd0H3A6sMLMN8fyxwAOSNgD7AH8U5/EcQeW+WdIPCIu7r5Xmrra/OHH5uwFYQPhFuBN4l5mtL1zjK0dfOXY+hq8cR123QxsyRx55pH3ve9+r6jtjxoy7xlCrpyyt1eq4q9T7RZgGXFISjI7jTB08z+uO2Rwxs+uB69tcm6+gpk2b9sLzfGWXrhDT5zNmzNiu3y677NJ3vNL4pVVkaQVbOl9aHXaxcnz22WdfeJ7bhtK+6Qqz9l5j3btNv9K9JqttkH/+Lt6D0qovbUuf559fepy3dS3MXDjuoHB0HGfnxIWjC0fHcTLMrMvd6inLpAnHXL0sbbQ0qdKzZs3arl+qVqfPoVmtHmRTp406XupXS0mtLqnOTepaqV+pb0n9bvvP1EYtHWteTW21cxzETabNZko+j/QzTD/bkgmmNH4X+MrRV46O4/TBhaMLR8dx+uDC0YWj4zgZtQ7eOztDY3MsOWk3ufLkdsXUHpnbEtO2ks0xtU3m7kDpdel8S/0GcfNpsk+WbIK1jsIll5/a8Uv2zVo73SC2w1p7Ya2tr4t7jTVOj5KrTd6W2hl/9rOfvfA8/660cdNqi2/I+MrRcZw++MrRhaPjOBmuVgeGRq2udaEpRdKkbSUXnSYVG9pF2bRVq2vV9pw2riWDqHlNbicl1bw2yqatWl1rIqg1MwziCtOFu1GqOqfP8+P0O1CKqe/CRayEC0dfOTqO0wcXji4cHcfpgwvHIRKObVKWDULTdYOoxE273PkY6Wsp7YaXkm/Umg9q34/aqIz8uFatHo8kF13vvNdGE41H9ElJrd62bVvfa0q72k3XdIGHDwaGRjg6jjM8+MrRhaPjOH1w4ejC0XGcPrhwHCLhWFueYJCSBE1jlFxmSi46TTbBvF9qZyxF+5TGb4oKgrK7UdP7U3IL2bp163ZtTXa6thEybVx+8nm0tSvWuiUNUIKAJtL3Ph+/KQomJ70u79d1tqcSXQpHSQsJdaSnAReZ2YqsfSZwOXAk8CRwmpk9Iul4YAWhJvU24GNmdku85jTgk3HMNWb28Xj+lYSiWnsRKhj+tpmNJPfajVDL5lozW1qa9/i+w47jTDl6GzI1j7GQNA24EDgRmAcsljQv67YE2GJmBwMXEGpUA2wGTjKzwwnVCa+IY84GzgcWmNl8YB9JC+I1XwAuN7NXA8uBz2b3+kPg72veBxeOjuOMohclM9ajgqOBjWb2sJltA64ETs76nAxcFp9fAyyQJDNbZ6EONcB6YFZcZR4EbLBQYhXgJmBRfD4PuDk+vzW9l6QjCRUJv1Uz8aGJkCm1NanVg1QYrK0cWFJdmlTpUnLeUqROKSFv2laqlZOP0aRuldTqmTNnbtfWxpWnNvHEIG4yTffOXWFqo33StpJrU9vXWatW559ROkbar/TdaWteqqVDtXo/4NHkeAQ4pqmPhcJ9TwOzCSvHHouAdWa2NZZwPVTSAXG8UwiqN8A9se9K4G3Ay+JKcwvwReB3CBVTx8RXjo7jjGKAleMcSWuTx5nZUP2keC55i30kzSeo2u+Lc9sCnAVcBXwHeATo/ap8FHiTpHXAmwhlo58FzgauN7NUUBcZmg0Zx3GGgwETT2y2ct3qEWD/5Hgu8FhDnxFJ04HdCZspSJoLXAu8x8weSuZ4HXBd7HMm8Fw8/xjw9nj+pcAiM3ta0huAX5N0NvBSYIakn5rZOU0Td+HoOM4oOlSr7wQOkXQgYRV3OvCurM9qwobLHcCpwC1mZpL2ANYAy8zs9vQCSXub2SZJexJWhe+M5+cAT5nZ88Ayws41Zvbu5NozgKNKghGGSDi2qfXcRWafUvhgyUWnZFesdcMpXVeyOaY2wpLNMX1/Sra+3GWkyeY4HmF1bVx0SvbC2lrPpX6l8Wttq/n4pWw76ftdCi+dSJtjV+GD0Ya4FLiB4HZziZmtl7QcWGtmq4GLgSuiLfEpggAFWAocDJwr6dx47gQz2wSslPSaeG65mW2Iz48FPivJgG8D728796ERjo7jDA9d+jma2fXA9dm5TyfPnwHe0ee684DzGsZc3HD+GsKOd2k+lwKXjjHtsTdkJF0iaZOkHybnXi7pRkkPxr97jjWO4zhTg9rNmJ09iqZm5Xgp8GWCB3uPc4CbzWyFpHPi8Se6mlRtfZlB6tA0qdy1UTD5cVvVubatpFanx7X1uUsqcT6PJjWyNjqk33ETJbW6ViUu1fFuastNCaUom9K9m8wO+fcv7VerLg9Sf6hrdnbBV8OYK0cz+zZx5yghddq8jOBn5DjOToKvHNvbHPcxs8cBzOxxSXs3dYzb7Lnvk+M4Q8zOLvhqGPcNGTNbBawCiDtIxOfb9WuTeKKkdtQmlBikDk1TW67atm1rUqVL/WrV6pz0y5+/ztqdytL7Xdp1bZpHSWWtTRZbUqtr67iU2krjl6KJ2iSNmEg1OqUXW/1ip22EzL9J2hcg/t3U3ZQcx5lsXK1uLxx7TpvEv9/oZjqO4wwDLhwr1GpJf0FwrJwjaQT4fUKOtaslLQF+RB8fJcdxpi47u+CrYUzh2ORsSWVmi1q6sDm2iXwpRankYzS54QzirpPaEmuvG6RIV5Ntq2S/ysdvyjTTNjFwOn4+RsmVp6mgVF5cqpRINj1u0y9vy69LX08piW1KbaLdyRRQLhw9QsZxnAzfkAm4cHQcZxS+chxS4VirKg6S7LbJxaVWHSy1lVTzPLql1pWn1l2nNMfa96rWjJH3K70H6euudWPJVyupmlpSq0tqb1NbaYy8rdY8kc4//16VGEZBNIxzmmiGUjg6jjO5uHB04eg4TsaLwU2nBheOjuOMwoXjEBXYqrWBtSmAlR+3qR0Nza42JdthbeadUtsgyXRrbatN1+TX1dpn89fS5LI0iM0xdal55pln+o4N29sISzbHtF8+37St9H6UXJHS0MLSGDltwgSnSrLbqYyvHB3HGYWvHF04Oo6T4TbHwNAIx9pIjFq1um3N6ZI626RG1mbXGWv8WrW91ixQ+37U1tvJ+7WprV1S70tqdfqac1ebdPxSxp60Xz5GSfWvNQWUalOXzBO1NZKmarJbSQsJdaSnAReZ2YqsfSYhmfaRwJPAaWb2iKTjCaHKM4BtwMfM7JZ4zWnAJ+OYa8zs4/H8KwlFtfYi5KH9bTMbkfRa4KvAboRKhX9kZleV5u11qx3HGUVXiSckTQMuBE4E5gGLJc3Lui0BtpjZwcAFhBrVAJuBk8zscEKCmyvimLOB84EFZjYf2EdSL5z5C8DlZvZqYDnw2Xj+PwnlXecDC4EvxeqGjbhwdBxnFB1m5Tka2GhmD5vZNuBKQiWBlLSywDXAAkkys3UW6lADrAdmxVXmQcAGM3sitt0ELIrP5wE3x+e39u5lZhvM7MH4/DFCmsW9ShMfyt3q2siXkurSNrqlpHI3JaAdJBltbeRLKUFFSQVsUrkHKUGbvt9tI4Ga5l+6V/7Plu7+1kbj5LvVaVuqSm/dunW7frVqdakWT3rvUsncQaKVJoOOY6v3Ax5NjkeAY5r6WCjl+jQwm7By7LEIWGdmW2MJ10MlHRDHO4WgegPcE/uuBN4GvEzSbDN7sjeQpKNj/4dKE/eVo+M4oxhg5ThH0trkkZdE6Sft8yVnsY+k+QRV+31xbluAs4CrgO8AjwA9g+9HgTdJWge8Cfhx0tZLzn0F8LtmVvwFGJoNGcdxhocBNmQ2m9lRhfYRYP/keC7wWEOfEUnTgd2JRf0kzQWuJdgLX1jpmdl1wHWxz5mETZaeyvz2eP6lwCIzezoe7wasAT5lZt8b64X5ytFxnFF0aHO8EzhE0oGSZgCnEyoJpKSVBU4FbjEzixsma4BlZnZ7ekGvqJ+kPYGzgYvi8RxJPbm2jLBzTbz3tYTNmr+smfjQrBzbRMjk9quSbavJZtU2gqXkalMqjtVFsts2rjy1dsXSdaV51NpWB6nF3FQELLdv1hbfSp+X3vvSHHOBkN4vnVca0ZPfr7am+s7gyhNtiEuBGwhuN5eY2XpJy4G1ZrYauBi4ItoSnyIIUIClwMHAuZLOjedOMLNNwEpJr4nnlpvZhvj8WOCzsZjft4H3x/PvBH4dmC3pjHjuDDO7u2nuQyMcHccZDrpOdmtm1wPXZ+c+nTx/hj6lVszsPOC8hjH7Vigws2sIO975+a8DXx9k3i4cHccZhUfIDKkrT6kGcttEtU1qatukEbV1pWsjR/LjUlKH2mS3bVTnUtt4JMCodWNJ34OSy09t4on8c0nV4EGSRjTVuSl9ZqUIsGFw5QEXjuArR8dx+uDC0YWj4zgZnngi4MLRcZxRuHAcIuFYGz5YGzpX66LT1l6YtpVC50ptpWw7bex5ed82WYpKbYOEWjbNv9aNBertkenOav6epjbH2jDG0vi5q9DMmTP73nsQe3htXfaJDDP0ZLdDJBwdxxkefOVYESEjaX9Jt0q6X9J6SR+M518u6UZJD8a/e47/dB3HGW9qo2N2dgFas3J8FviImf2TpJcBd0m6ETgDuNnMVkg6BzgH+ETtjQdx5WlSh0pqadua0G3ccPJ+qapVq5rn82+TeSfvW5vQtm2ETDpGbULeQdTq2ho4qQpYGqMU6ZKOkboGQdlFpwtTyM4cITOVGXPlaGaPm9k/xef/AdxPSDGU5mC7jJA2yHGcnQBfOQ5oc4z5044Avg/sY2aPQxCgvUDwPtecCeRpjBzHGWJ2dsFXQ7VwjOl//gr4PTP799olvpmtAlbFMfwdd5whp+vY6qlKlXCUtAtBMP65mf11PP1vkvaNq8Z9CWnHW5MK21o72iDFsdq44ZQy6tSOkdofxxq/yZ5acguptbHVuuuU2mozHQ0yRhubY8leWCLtN0g4ZW0WpJKNtI3L0mSGEvrKsW63WoSUQveb2R8nTWkOtvcC3+h+eo7jTAZuc6xbOb4R+B3gB5J6uc/+F6Fk4tWSlgA/ok/KIcdxpiY7u+CrYUzhaGb/QP8aDwALGs6PySCuPE1uJ7Wqc35ccuVpE/mSq87pcamtNitPyY2lpL6Voi26qPHdJrlw2wiZUpGrpvHy49qEym3b2rrh1GblmUiV24WjR8g4jpPhGzIBF46O44zCV45DlOy2TW2YQXaJ0+NZs2Y19mujVufqcUmtTu9dip4pqb0ptSpgrboJze/3IMlumzwHSvOopfSPW9rJTp/n/dLj0hi1QmOQaJ+m96DNe9MVLhy9+qDjOH3ocrda0kJJD0jaGEON8/aZkq6K7d+PwSZIOl7SXZJ+EP8el1xzmqR7Y76HzyfnXynp5th2Wyzt2mt7b8wF8aCk9zIGLhwdx9mOLhNPSJoGXAicCMwDFkual3VbAmwxs4OBC4DPxfObgZPM7HCCu+AVcczZwPnAAjObD+wjqbc5/AVC+dVXA8uBz8ZrXg78PnAMcDTw+2Mly3Hh6DjOKDpcOR4NbDSzh81sG3AlIS9DSpqn4RpggSSZ2TozeyyeXw/MkjQTOAjYYGZPxLabgEXx+Tzg5vj81uRevwncaGZPmdkW4EZgYWniQ7MhU1tgq1S/OD1ObXv5ccnm2CZjT8m+ueuuuzbOo5SVp+S+M4A609jWtStPbe3uQepWNyWZzXdSS3Wrm9pKhbjytvS6PGNP7a5urTtQqRDXREbMDLBbPUfS2uR4VQwZ7rEf8GhyPEJYvdGvT6xz/TQwm7By7LEIWGdmW2N960Oj+j1CSHrT+8LdE/uuBN4GvCyuNPvNY7/SCxsa4eg4znAwYPTLZjM7qtDeT6Lngxf7SJpPULVPiPPbIuks4CrgeeC7hNUkwEeBL0s6A/g28GNC2sWaeWyHC0fHcUbR4W71CLB/cjwXeKyhz4ik6cDuwFMAcUPlWuA9ZvZQMr/rgOtinzOB5+L5x4C3x/MvBRaZ2dOSRoBjs3ncVpr40LjytKkhk6t5Xbjh1KrVtcluSyp3aY4lN44m95QSg7zftUkj2iQB6cJckL/mUqLaJlW6pH7nY6THtWp1yVVtPKJ4uqZD4XgncIikAwmruNOBd2V9enka7gBOBW4xM5O0B7AGWGZmt6cXSNrbzDbFTZWzgXfG83OAp8zseWAZcEm85AbgfyebMCfE9kZ8Q8ZxnFF0tSFjZs8CSwnC6X7gajNbL2m5pLfGbhcDs6Mt8cOEqgLE6w4GzpV0d3z08saulHQfcDuwwsw2xPPHAg9I2gDsA/xRnMdTwB8ShPWdwPJ4rhFXqx3H2Y6uwwfN7Hrg+uzcp5Pnz9AncY2ZnQec1zDm4obz1xB2vPu1XcLPV5Jj4sLRcZxReITMkArHkitPmzDDUtsgSU+brsv71Wb9KbkRlWyOqd0r/4Vv4+ZTmylnkMw+tYW+UkqrlVJ4X8km2PRelWyHXbynbW3qk5ngNsWF45AKR8dxJhcXji4cHcfpgwvHIRKOtapF10lJB0m+Whs5UuviUlvrufRFzV9nk/rZNrlwm2S0/cZsQ9NrKbny1LYNMkZJpW+TBamt61RtUtwd5cVQAqGGoRGOjuMMD57s1oWj4zh98JXjFBSOJRWnNglqrZraNpFsm7Kn0JwUtpR8NZ9jU2LW2homed/a96NE7Xs/Hipx0271IKpzrXmi1pOi1kNiMhNPuHCcgsLRcZzxxW2OAReOjuOMwoWjC0fHcfrgwnEKCsdSlENt9pTJjIYouWPUZmdJ7VL566yt79w030HaSu41+bya+pWuafo820a3tP3ca+21bYuRNdma834TGUnju9VTUDg6jjO+uM0xMGbKMkmzJP2jpHtipa8/iOcPVKgU9qBC5bAZY43lOM7UoMvqg1OVmpXjVuA4M/uppF2Af5D0t4S8axeY2ZWS/oRQQeyrXUyq5EqRPs/VqzSZ6bZt27Zra4o+KUUh5MloS3PsmlRtKkVldDF+iVoXl1Jdl/Rzys0MtW44pRoybV10mvqVKLlt1SZArk3EPIhJpmt2dsFXw5grRwv8NB7uEh8GHMfP86ZdRihy4zjOToCvHCttjrH27F2ErLwXAg8BP4lZfqFQySvWdzhzx6fqOM5E0HWy26lKlXA0s+eA18aaDtcCh/Xr1nDtKmAVgKSd+6fGcXYSdvZVYQ0D7Vab2U8k3Qa8HthD0vS4euxXUWwgasPKal158gJKtW4hJZtmelxrixuPEMfJ+uIOEsaYkraVbILj4crTFE45CLWuWW0zNTW58pTuNd50+R2TtJBQR3oacJGZrcjaZwKXA0cCTwKnmdkjko4HVhBqUm8DPmZmt8RrTgM+GcdcY2Yfj+dfQTDz7RHbzjGz6+N+yUXA6why73Iz+2xp3jW71XvFFSOSfgF4M6FQzq2ESmEQKod9Y6yxHMeZGnRlc4wmuQuBE4F5wGJJ87JuS4AtZnYwcAGhRjXAZuAkMzucIGOuiGPOBs4HFpjZfGAfSQviNZ8iFPE6glDp8Cvx/DuAmXGsI4H3STqgNPean6J9gVsl3Uuo2nWjmX0T+ATw4VgxbDahgpjjODsBHW7IHA1sNLOHzWwbcCVwctbnZMJqD8Im7wJJMrN1FupQA6wHZsVV5kHABjN7IrbdBCzqTR3YLT7fnZ9rtAbsGuti/wJhJfrvpYmPqVab2b3AEX3OP0x44Z1Tq1KWXD9Kbj6pyl1Sv9tGW3TRVqKtSl9LkxvRIK48TePl1KrVXSSqLdE2iXKbejsl97HaejtDlOx2jqS1yfGquM/QYz/g0eR4BDgmG+OFPmb2rKSnCQuuzUmfRcA6M9saF2SHxpXfCMFTpucD9RngW5I+AOxK0HQhCN2TgceBXwQ+ZF6a1XGcQRlgt3qzmR1VaO8nxXPJW+wjaT5B1T4BwMy2SDoLuAp4HvguYTUJsBi41My+KOkNwBWSXkVYyD0H/BKwJ/AdSTfFRV5fJs7C6zjOlKFDtXoE2D857rd5+0KfqPbuDjwVj+cSPGTeY2YPJfO7zsyOMbM3AA8AD8amJcDVsc8dwCxgDvAu4O/M7Gdmtgm4HSgJ9eFZOdYG1ZfUvNrSm7WlPGvbctW8pLanx7l6lb7uvK1pHm3VyFpKO9K1n0Xb5BVtkt12YWbooqzqZCaq7YIOv0d3AodIOhD4MWGT5F1Zn9WEDZc7CJu8t5iZxY3gNcAyM7s9vUDS3ma2SdKewNnAO2PTj4AFwKWSDiMIxyfi+eMkfZ2gVr8e+FJp4r5ydBxnO2pXjTUCNLr6LQVuIHi5XG1m6yUtl/TW2O1iYHa0JX4YOCeeX0oIPDlX0t3xsXdsWynpPsIKcIWZbYjnPwL8D0n3AH8BnGFhohcCLwV+SBDYfxb3UxoZmpWj4zjDQ5caiJldD1yfnft08vwZgqtNft15wHkNYy5uOH8f8MY+53/a7x4lXDg6jjOKyQo0GCYmTTgO4i7RRMnOVbJLleyKTdlk8rbafqkLUX5cctWoLaLVNga21j0lpSmBbT9qozlq7ZYT6b5Um4S4xCBzbGqbzPhmj632laPjOBkD+jnutLhwdBxnFC4ch0itrg24b+sW0qQGt01eUas650l3Sy466ZzTZAQlF5oStckgamuktJ1HaU7pGKV5dFFzukSp5nQbN7NBordKrllN44+38HLh6CtHx3H64MLRhaPjOBnmyW4BF46O4/TBV46TKBxzu04pG0lTMtCcks2nKfvLIOGDqS0xLaaU2xxTO2OtjRG2t2OWXnPJDtgmnK2L0LZSWF3pM6uts13r8lMbTjnI55J+FrkdOp1/el3eL/2O5N+XJhex8Q4NLeHC0VeOjuP0wYWjC0fHcfrgwnGIXHlKGWm6yNjT5C5RGwWTH6eqUV4TJFWra80ApXnl45fcnprU2bZRH7X/JLXuNW3vVVs7aLxryNTOqzbZcn5cW597PIWXO4EHfOXoOM4ofLfahaPjOH3wleMUF461ySVKbYNEyNSq1bVqb07Tjnpp97RUj6REeq/8mlJbU7/a5BWlMdqWZm3bllJSWWt3w2tL/taqy5O5enPhOMWFo+M43eM2x4ALR8dxRuHC0YWj4zh9cOE4pMKx1g2nlPmkNitK2+JYpaw8tRlecrqIVGlyocltgKX3sSnqo1SLOSd9LaXEvaXPrCmqpJQFKW+rdX+pzcZU60ZUsnlPBZtjl/eWtBBYCUwDLjKzFVn7TOBy4EjgSeA0M3tE0vHACkJN6m3Ax8zslnjNacAn45hrzOzj8fwrgMuAPWLbObFMA5JeDfwpsBuhpOt/iyUa+uIFthzH2Y4uC2xJmkYobnUiMA9YLGle1m0JsMXMDgYuINSoBtgMnGRmhxOqE14Rx5wNnA8sMLP5wD6SFsRrPkUo4nUEodLhV+I104GvA/8zXnMssP0vaUa1cJQ0TdI6Sd+MxwdK+r6kByVdJWlG7ViO4ww3XQlH4Ghgo5k9bGbbgCuBk7M+JxNWewDXAAskyczWmVmvxvV6YFZcZR4EbDCzJ2LbTcCi3tQJK0MI9a97158A3Gtm98TX96SZFet+DKJWf5BQWrF3488BF5jZlZL+hCD9v1o7WCkpaa27RK7W1UZRNCV4yNtqXXlK7jS1dY5LtI1uaVJtx7qu6Us/SC2bWneg0mfWlDQ4TyBcij6prbeTtpXU9pLKXesOVJtgo9Q2RMlu50hamxyvMrNVyfF+wKPJ8QhwTDbGC33M7FlJTwOzCSvHHouAdWa2NZZwPVTSAXG8UwiqN8BngG9J+gCwK/DmeP6XAZN0A7AXcKWZfb70wqpWjpLmAu7N31kAABXPSURBVP8duCgeCziOIOUhSP1TasZyHGf4GWDluNnMjkoeq7Kh+v0K55K32EfSfMJi7H1xbluAs4CrgO8AjwC9X63FwKVmNhd4C3CFpJcQFoK/Crw7/n1boor3pVat/hLwcYIRE4JU/4mFgt0QpPd+/S6UdKaktdmvi+M4Q0ov2W3No4IRYP/keC4/V3VH9Ym2wd2Bp+LxXOBa4D1m9lAyx+vM7BgzewPwAPBgbFoCXB373AHMAubEe/y9mW02s/8k1NF+XWniYwpHSb8FbDKzu9LTfbr2XYeb2arer8pY93IcZzjo0OZ4J3BI3KOYQdgkWZ31WU3YcAE4FbjFzEzSHsAaYJmZ3Z5eIGnv+HdP4GyiVgv8CFgQ2w4jCMcngBuAV0v6xSiA3wTcV5p4jc3xjcBbJb0l3mg3wkpyD0nT4+qx369BkZKrTa0LTcm1JLcbpW2l5KWltiY7Y36vUqLX2kJi6b1LiYG7sF/VtpXccPLXks6/ZN+sdeXZunVr3/P5ca2rTckmmH/uqY0znUfp3l3UE59MurJpRhviUoJwmgZcYmbrJS0H1prZauBigvq7kbBiPD1evhQ4GDhX0rnx3AlmtglYKek18dxyM9sQn38E+JqkDxEWbGdYeDFbJP0xQVgbcL2ZrSnNfUzhaGbLgGUAko4FPmpm75b0lwQpfyVB6n9jrLEcx5kadLnhE/0Mr8/OfTp5/gzwjj7XnQec1zDm4obz9xEWdP3avk5w56liR/wcPwF8OEr72QTp7zjOTkCHavWUZaAIGTO7DbgtPn+Y4MPUilq3DWiu4VxSq3NVtI0rT8lFp1TzJj0uRc/U1u7O59jFl7LkOtUmqWopeXHTffN71SYe7kKtHiQbU+m7mR6Xsv7UmlqGQa1+MQi+GoYyfNBxnMnFk926cHQcpw++cpxE4VgqS9o24UNJrW7ahc7vVRojPa7dNR8keqZpN3yQRAhN0Si5upZ++WujPnJq/4FKkR1t1OpB6v60SZhbUqtLO9mlOU5U/ZeumApzHG985eg4zna4zTHgwtFxnFG4cHTh6DhOH3xDZohsjm1sT4PY4ppcdAaxOda68qQ2qtq60qXxS/bTtnWrSxEhTe9jKZNSbdt4JyiutSW2LaxW2zbId7OLjD1d4mp1wFeOjuOMwoWjC0fHcfrgwnEShWNbl45S8tJSLelaFa2UeKIp2cQgySXaJKXI+zVFDEFzktmSWp27RNXW7GmT3HUQc0rtZ1arEteaEroYf5CEvE3qfheRS21x4egrR8dx+uDC0YWj4zgZvWS3L3ZcODqOMwpfOQ6RcCzZpdK2WpeIki2xjcsPNNsZ29ocS5lsStlZmupKQ3MC2i5ceUqfS8k+VrIrNn22pXmUPrO2rjyltvRzL4W9luy4aZLcUsLc0jyGtMDWTsvQCEfHcYYHF44uHB3HyXAn8MBQuvLkbW0yq9TWqMlV59pkt6UaKaU6zbXUqlAlFbCkVpey8jSpxKX3tPR5ltTqNq48tf2gWU1t665Te+9crX7mmWf6Poft1eymLD/5+FNJrZa0EFhJqCFzkZmtyNpnApcDRwJPAqeZ2SOSjgdWEGpSbwM+Zma3xGtOAz4Zx1xjZh+P519BKBW9R2w7J5ZpIGm/D/iMmX2hNO8d/y92HGeno6vSrJKmARcCJwLzgMWS5mXdlgBbzOxg4AJCjWqAzcBJZnY4oU7VFXHM2cD5wAIzmw/sk9Sg/hRwtZkdQSjU9ZXsXhcAf1vzHrhwdBxnFB3WkDka2GhmD5vZNkJBvpOzPicTVnsA1wALJMnM1plZr6rpemBWXGUeBGwwsydi203Aot7UCRVSIdS/fqEqqqRTgIfjWGMyNGp17W51rZo3Hipak+pVSjyRU9qFbpOsoZTwt7RrXvrVbxPd0na3uk3iifw1t6kN0zbxRMmMUUrEnKrO//Vf/7VdW3pcKkHbRenXGjq2Oe4HPJocjwDHNPWJpVyfJhTt25z0WQSsM7OtsajfoZIOiOOdQlC9AT4DfEvSB4BdgTcDSNqVUBTweOCjNRP3laPjOKMYYOU4R9La5HFmNlS/1UAueYt9JM0nqNrvi3PbApwFXAV8B3gE6P0yLQYuNbO5wFsI9bBfAvwBcIGZ/bT2PfDdasdxRjHAynGzmR1VaB8B9k+O55KoulmfEUnTCerwUwCS5gLXAu8xs4eS+V0HXBf7nAn0ltVLgIWxzx2SZgFzCKvVUyV9nrBZ87ykZ8zsy00Td+HoOM4oOlTb7wQOkXQg8GPCJsm7sj6rCRsudwCnAreYmUnaA1gDLDOz29MLJO1tZpsk7QmcDbwzNv0IWABcKukwYBbwhJn9WnLtZ4CflgQjDJFwrC3CVOsWUoqQaVu/uMmGl9uGamsPt7XTlepu19ocS7Sxfbb9zMbbTtwU3TJIhExtstva+tYlV570ecnFajxdebq0OUYb4lLgBoJrzSVmtl7ScmCtma0GLiaovxsJK8bT4+VLgYOBcyWdG8+dYGabgJWSXhPPLTezDfH5R4CvSfoQQTU/w1q+mKERjo7jDA9dCt/oZ3h9du7TyfNngHf0ue484LyGMRc3nL8PeOMY8/nMmJOmUjhKegT4D4Je/6yZHSXp5QSD6AEEg+g7o6HUcZwpjkfIDLZy/A0zS7fWzwFuNrMVks6Jx5/oYlJdqNX5GE2q6CDuKU0qd65GpypU7WvJx69NyFuqrV2bAKNEG3MHNLudtI1u6UIlHg+1usmVp61aXUqKO1UjZKYqO+LKkzpuXkbwNXIcZyegQyfwKUvtytEIjpUG/KmZrQL2MbPHAczscUl797swbrPnvk+O4wwpnuw2UCsc32hmj0UBeKOkf669QRSkqwCicHUcZ8jZ2VeFNVQJx158Y/QrupYQL/lvkvaNq8Z9gU2D3LiL8MGSDSy3G9Umuy21tXWNaZpjPv/UXpjeu1S3um1t7drEuimDvN9NbV1k1OnCbtm25nmtK08pfLC2+NZkufJMxPhTgTH/wyXtKullvefACcAP+bnjJvHvN8Zrko7jTCxuc6xbOe4DXBtXF9OB/2NmfyfpTuBqSUsIXumj/JQcx5l6vBgEXw1jCkczexh4TZ/zTxLCdDqh1mUkbWub4aU2m0xJNU/75WppqiYNooo2qcG17jp5W6lfGzefQV5LrUrcJpHsVFCr836pKl2rVpdep6vV449HyDiOMwrfrXbh6DhOH3zl6MLRcZwMtzkGhjITeK2bT96vtqZ1bS3mnKYxSq4wg2TsqQ3968Lm2NSvX98epfe7VM+51mbXtjhWm0zjbcMHa+2RuV2xlCW8qaiW162eXHzl6DjOKFw4unB0HKcPviEzRMKxFCFTG23RFGGSH9dGdtQySLaaEqnKXetqU+vmkyfFTY9LCXPTe5VeZ64qNkWLlNTqWvea8VCr2xRWK7WVsveUCmc1uR7lx1Ml2e1UZmiEo+M4w4MLRxeOjuP0wYXjFNmtblJTa5PRQnPiidokF/lxbf3p0hglmlTsvK028URJrZ4xY8Z2bfmYTZQiTpoSuOa7uG3qUXcRGVXadS7typfMAqX3o82O/SDJnLumy/ElLQRWEmrIXGRmK7L2mcDlwJHAk8BpZvaIpOOBFYSa1NuAj5nZLfGa04BPxjHXmNnH4/lXEPLL7hHbzjGz60tjNeF1qx3HGUVXiSckTQMuBE4E5gGLJc3Lui0BtpjZwcAFhBrVAJuBk8zscEJymyvimLOB84EFZjYf2EdSL5T5U8DVZnYEoVDXV0pjlXDh6DjOdvSS3dY8Kjga2GhmD5vZNuBKQhWBlLSqwDXAAkkys3W9dInAemBWXGUeBGwwsydi203Aot70gd3i892JNbILYzXiNkfHcUbRoVq9H/BocjwCHNPUx0Ip16eB2YTVXo9FwDoz2xpLuB4q6YA43ikEdRngM4SqBR8AdgXe3GdOL4xVmvjQCMc2NsdSMtrcJpjadUoZb1JbXMmumN47twmW5tgmyiGfR8nNp+m15a8ztTOWXKJq34NS1EdTBEh+XOvKUxv9lB+XMjq1cQcqzavW5Sc/rrWHD5HNcY6ktcnxqpj9v0e/L1A+eLGPpPkEVfuEOLctks4iVD99HvguYTUJsBi41My+KOkNhHrYrzKz5/uNVWJohKPjOMPDAMJxs5kdVWgfAfZPjucSVd0+fUYkTSeow08BSJoLXAu8x8weSuZ3HXBd7HMmoWw0BPvlwtjnDkmzgDnApqaxmnCbo+M421G7GVMpQO8EDpF0oKQZhE2S1VmftKrAqcAtZmaS9gDWAMvM7Pb0gl5BP0l7AmcDF8WmHxHzzEo6DJgFPFEaq4mhdOWpjZAZpBZzrvY1jVFSzdM5llxt2tbWbjJwDzLH2giZkntN2jcdY5AIma7V6jaqc6mtZO4ofWalBBtt3czaJGIeIrV6rHGelbQUuIHgWnOJma2XtBxYa2argYsJ6u9Gworx9Hj5UuBg4FxJ58ZzJ5jZJmClpF4S7uVmtiE+/wjwNUkfIqjmZ0RBWxqrL5pIZ8+0+mDJ1lcb6rbLLrtUjzFz5s83ptLrZs2atV2/9Dhva7JHTkXhmL4HpffRhePECsda22RpDDOrq5jWwIwZM2yvvfaq6vvYY4/dNYZaPWVxm6PjOKPwCBkXjo7jZHjiicDQCMc2NseSjbFU9CqlVH86/4I0qZul60oZZNqmhSqp1elxqb51Sa2udeUphculqnTbrDxNqu4gtrgu1Ora72atWl3rDjRIIuauceE4RMLRcZzhwYWjC0fHcfrgyW6niCtPrua0Gb82yqakKtaqmyklVbHtr3N679romVJS3FytLpkMUkrmg9oaMl1EtzSNlx/XRmHVfnfyvrVuOLVmI092O7lU/RdI2kPSNZL+WdL9kt4g6eWSbpT0YPy753hP1nGciaFDJ/ApS22EzErg78zsUOA1wP3AOcDNZnYIcHM8dhxnJ8CFY4VaLWk34NeBMwBi2qFtkk4Gjo3dLgNuAz5Re+OSWl3rHN1FlE3uHJ2qfbkqWkoyWzuPpn79jvvdd6y2pp3svF9p572NWl2byGEQB+um93GQz73pezXI96/N+G2TKLe513iwswu+Gmr+Cw4CngD+TNI6SRdJ2hXYx8weB4h/9x7HeTqOM4H4yrFOOE4HXgd81UJ23f/HACq0pDMlrc3SGjmOM6RYt8lupyw1wnEEGDGz78fjawjC8t8k7QsQ//YN4DazVWZ21M4af+k4OyO+cqywOZrZv0p6VNKvmNkDhHRA98XHewlFa94LfGNHJlJrm2tjV4RmF5fcvtaUuAHKNryauQ/SVqLWBlnr8lOKsmlbSKzWxaWte03tGLV2ui7skaXzbV5n7b3Gg51d8NVQ6+f4AeDPYz62h4HfJaw6r5a0hJBD7R3jM0XHcSYaF46VwtHM7gb6qcUL+pxzHGcK82JQmWsYmgiZlJKKXaodXVIHa9XN2rYuGA+1urZfGxNBiS5U1kGu29F5DNMYtap5m3m0xYWjx1Y7jtOHnX0nugYXjo7jjMJXjl5gy3GcjFo3nloBKmmhpAckbZQ0ykda0kxJV8X278d61Eg6XtJdkn4Q/x6XXHOapHslrZf0+eT8KyTdGgNW7pX0lqRtWbzHA5J+s9M3YkcfhII3/vCHP8bx0cX/6bRp06oehCJZpbGmAQ8RIu1mAPcA87I+ZwN/Ep+fDlwVnx8B/FJ8/irgx/H5bIKHzF7x+DJgQXy+CjgrPp8HPJI8vweYCRwY5zStNHdfOTqOM4oOV45HAxvN7GELeRmuBE7O+pxMEHAQgkwWSJKZrTOzXo3r9cAsSTMJgnaDmT0R224CFvWmDuwWn+/Oz2tknwxcaWZbzexfgI1xbo24zdFxnFF0uCGzH/BocjwCHNPUx0Ip16cJq8PNSZ9FwDoz2xpLuB4a1e8R4BTCqhTgM8C3JH0A2BV4c3KP72Xz2K808YkWjpuB/wvMYfsXPhkMwxzA55Hj89ieQefxyg7ueUO8bw2zsrwJq8xsVXLcz08sX3IW+0iaD3wOOAHAzLZIOgu4Cnge+C5hNQmwGLjUzL4o6Q2EetivqpzHdkyocDSzvQAkrbVJjrUehjn4PHwewzgPM1vY4XAjwP7J8Vx+rurmfUYkTSeow08BSJoLXAu8x8weSuZ4HXBd7HMm0IsXXgIsjH3ukDSLIOhr5rEdbnN0HGc8uRM4RNKBMfz4dGB11mc1IT8DwKnALWZmkvYA1gDLzOz29AJJe8e/exI2dC6KTT8iRu5JOgyYRUi5uBo4Pe6MHwgcAvxjaeJuc3QcZ9yINsSlBFV9GnCJma2XtJyw070auJig/m4krBhPj5cvBQ4GzpV0bjx3gpltAlZKek08t9zMNsTnHwG+JulDBLX5DAs7R+slXU1ImPMs8H4zKxan0mQ4e0o6M7NLvCjn4PPweUyVebwYmRTh6DiOM+y4zdFxHKcPEyocxwojGsf7XiJpk6QfJucmvLSspP1jaNP9Mezpg5MxF0mzJP2jpHviPP4gnj8whm89GMO5Zow1VkfzmRbDvb45WfOQ9EgMU7u755oySd8RL4M8JEyYcJQ0DbgQOJEQyrNY0rwJuv2lxO39hMkoLfss8BEzOwx4PfD++B5M9Fy2AseZ2WuA1wILJb2e4Et2QZzHFoJbxETwQUK53x6TNY/fMLPXJq4zk/Ed8TLIw8J4xFA3hBm9AbghOV5G2KKfqPsfAPwwOX4A2Dc+3xd4YKLmkszhG8DxkzkX4BeBfyJELWwGpvf7vMbx/nMJ//DHAd8kOOtOxjweAeZk5yb0cyGEvf0LcS9gsubhj/CYSLW6XxhRMXxnnJnU0rIx9OkI4PuTMZeoyt5NKIx2IyEQ/ydm1ivcPVGfz5eAjxMiHSCEjU3GPIwQdnZXdCqGif9cvAzyEDGRwnHg8J2dFUkvBf4K+D0z+/fJmIOZPWdmryWs3I4GDuvXbTznIOm3gE1mdld6eqLnEXmjmb2OYPZ5v6Rfn4B75uxQGWSnWyZSOA4cvjPOVJWW7RpJuxAE45+b2V9P5lwAzOwnwG0EG+geMXwLJubzeSPwVkmPELK1HEdYSU70PLCY/cWCg/G1hB+Mif5cdqgMstMtEykca8KIJpI0ZOm97GBp2RokiRANcL+Z/fFkzUXSXjE0C0m/QMhccj9wKyF8a0LmYWbLzGyumR1A+D7cYmbvnuh5SNpV0st6zwkJDn7IBH8uZvavwKOSfiWe6pVBnvDvqsPEbchYMCa/BdhAsG99cgLv+xfA48DPCL/OSwi2rZuBB+Pfl0/APH6VoCLeC9wdH2+Z6LkArwbWxXn8EPh0PH8QId50I/CXwMwJ/IyOBb45GfOI97snPtb3vpuT9B15LbA2fjZ/A+w5GfPwh3mEjOM4Tj88QsZxHKcPLhwdx3H64MLRcRynDy4cHcdx+uDC0XEcpw8uHB3HcfrgwtFxHKcPLhwdx3H68P8BYsLpIsgCOJoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "im = plt.imshow(U, cmap='gray')\n",
    "plt.colorbar(im)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 236,
   "metadata": {},
   "outputs": [],
   "source": [
    "laplace = torch.tensor([[0,  1,  0],\n",
    "                        [1, -4,  1],\n",
    "                        [0,  1,  0]], dtype=torch.float32).unsqueeze(0).unsqueeze(0)\n",
    "\n",
    "UU = torch.tensor(U, dtype=torch.float32).unsqueeze(0).unsqueeze(0)\n",
    "conv_U = nn.functional.conv2d(UU, laplace)/(h**2)\n",
    "pad = nn.ZeroPad2d(1)\n",
    "conv_U = pad(conv_U)\n",
    "conv_U = conv_U.numpy()[0,0,:,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 237,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'conv u')"
      ]
     },
     "execution_count": 237,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA24AAAEoCAYAAAAgz6I/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9f9g1VX3e+7kFXmg0CoIggogWUkWjMbxFPZ5EK0JIW4sJKmCOQosh0WDbJP6AS7GWkiNqjOUcbdo3SEBsI5YTK0QSRAixVWN5uVATMAGkKq8Q4QWkUcPv7/lj5oF5N3vtZ82emb1nr+f+XNdzPTNr1lqzZvZ+7mfWrHt9lyICY4wxxhhjjDHj5XHLboAxxhhjjDHGmNm442aMMcYYY4wxI8cdN2OMMcYYY4wZOe64GWOMMcYYY8zIccfNGGOMMcYYY0aOO27GGGOMMcYYM3LccTPGGGPMaJF0lKS/lnSTpFOnHN9V0oX18a9IOnDxrTTGmOFxx80YM/eDkaQjJF0j6S/q369olDlW0tclXSfpA430Z0i6oj52laT9G8cOkPQ5Sd+QdL0fwIzZ2EjaCfgo8PPAIcDxkg6ZyHYScHdEHAR8GHj/YltpjNmILOOlkjtuxmxwOj4YbQdeFRE/CZwAXFDXuSfwQeDwiHgusI+kw+syvw18PCKeD5wBvK9xno8DH4yI5wCHAbf3erHGmFXjMOCmiLg5Iu4HPgkcPZHnaOD8evsi4HBJWmAbjTEbjGW9VHLHzRgz94NRRFwbEbfW6dcBu0naFXgWcENE3FEf+zxwTL19CHBFvf2na+eqBW/niLgcICJ+EBE/6vNCjTErx37ALY39bXXa1DwR8SBwD7DnQlpnjNmoLOWl0s5dChtjFstRRx0V27dvb13ummuuuQ64t5G0JSK21NvTHoxeNFHFDg9GktYejJqNOQa4NiLuk3QT8OzaFrANeDWwqc73tTrv2cAvAD9ej9D9BPB9SX8IPJOqs3dqRDzU+oKNMQtlIG0CmPaQExP7OXmMMRuUefTpmmuuuSwijpqRpa9np1a442bMCrF9+3a2bt3aupykeyNic+rwlLRWD0aSnktlATgSICLulvRm4ELgYeBLVKNwAG8DPiLpROALwHeBB6n06GeAFwLfqcueCHxs9tUZY5bNQNoE1cPQ0xv7+wO3JvJsk7Qz8CTgrtaNMcYUyTz6JOnZkpqFRvFSyR03Y1aMiN5fJHd6MKqDi3waeGNEfLPRzkuAS+o8JwMP1em3Ar9Ypz8BOCYi7pG0jWrE7ub62H8DXow7bsasBANoE8DVwMGSnkn1kuc44PUTeS6mmmP7ZeA1wJUxUGOMMavJHJKwfYwvlTzHzZgVIyJa/6zDIw9GkjZRPRhdPJFn7cEIGg9GknYHPgucFhFfbBaQtHf9ew/gLcA59f5ekta05zTg3EY79pD0lHr/FcD12TfGGLNUBtCmtTlrpwCXAd8APhUR10k6Q9I/q7N9DNiztmj/BvCY6G7GmI1N39pEh2enLtfhETdjVoy+XyTXvuu1B6OdgHPXHoyArRFxMdWD0QX1g9FdVAIF1QPVQcDpkk6v046MiNuBsyW9oE47IyJuqLdfDrxPUlBZJX+tbsdDkt4GXFFP3r0G+L1eL9YYMxhDDXJFxKXApRNp72ls3wu8dpCTG2OKYGTPTnMjuwmMWR0OPfTQ+MpXvtK63C677HLNOkP+xhgzN9YmY8xYmUefxqpNHnEzZsXwyxZjzBixNhljxkop+uSOmzErRiniY4wpC2uTMWaslKJP7rgZs2KUIj7GmLKwNhljxkop+uSOmzErRiniY4wpC2uTMWaslKJP7rgZs0K0CFNrjDELw9pkjBkrJemTO27GrBiliI8xpiysTcaYsVKKPrnjZsyKUYr4GGPKwtpkjBkrpeiTO27GrBiliI8xpiysTcaYsVKKPrnjZsyKUYr4GGPKwtpkjBkrpeiTO27GrBAlTbA1xpSDtckYM1ZK0id33IxZMUoRH2NMWVibjDFjpRR9csfNmBWjFPExxpSFtckYM1ZK0afHLbsBxhhjjDHGGGNm4xE3Y1aMUt4aGWPKwtpkjBkrpeiTO27GrBAlTbA1xpSDtckYM1ZK0id33IxZMUoRH2NMWVibjDFjpRR98hy3dZD0DyRdK+lvJf1LSf9R0umZZbPzztm2kHTQtHNJerOk70n6gaQ9Jb1U0o31/quHapMZnrU3R21+jCmZSb1bdns2KtYmY8xYKUWbPOK2Pu8AroqIF7YtGBG/urYt6eXAJyJi/x7bljrXLsDvAC+OiK/VaWcAH4mIs4c4v1kcYxYUYyR9C3hTRHx+Qed7jN6Z5WBtMsaMlVL0yR239XkG8MllN6Il+wC7Adc10p4xsW9WlFLEx2xMJO0cEQ/2WOU0vTNLwNpkjBkrpeiTrZIzkHQl8I+Aj9QWnJ+QdJ6kM+vjL5e0TdJvSrpd0m2S/nmj/HmSzpT0eOCPgafV9fxA0tMkPU7SqZK+KelOSZ+S9OQZ7Xl7fY5bJf2LiWNr5/oJ4K/r5O9LulLSN4FnAZfU595V0pMkfayu77t12Z3quk6U9D8k/bakuyX9L0k/3zhXl7JPlvT79TXcLem/NY79U0lflfR9SV+S9PzEfTiwtonu3Ei7StKb1vtMV515rEiliJUZP5IuAA7gUa15R+Pv9SRJ3wGuXNPOibLfkvTKejtLG6fp3cCXaBJYm8wYkfR0SX8o6Y5aSz5Spz9O0rslfbt+fvu4pCfVx9Y06wRJ35G0XdK76mNPk/R3TT2S9MI6zy5Tzv/IM2O9/xjtM8NTkja54zaDiHgF8N+BUyLiCRFxw5RsTwWeBOwHnAR8VNIeE/X8EPh54Na6nidExK3AvwReDbwMeBpwN/DRaW2RdBTwNuAI4GDglYk23wA8t97dPSJeERF/H/gO8Kr63PcB5wMPAgcBLwSOBJodnxdRPRDtBXwA+Jgk1ce6lL0A+LG6jXsDH66v76eBc4FfAfYE/hNwsaRdp13nRqYkATJlERFvYEet+UDj8MuA5wA/l1FVljZO07sOzTcdsTaZMVG/UP4j4NvAgVTPaWsOqhPrn39E9WL7CcBHJqr4P4F/ABwOvEfSc+pnty8DxzTyvR64KCIeGOI6TD+Uok3uuHXnAeCMiHggIi4FfkD1h57DrwDviohtdWfqvcBrmiNJDV4H/H5E/GXdEXzvvA2WtA9VR/JfR8QPI+J2qg7UcY1s346I34uIh6g6avsC+3Qsu29d9lcj4u76nv1ZXeaXgf8UEV+JiIci4nzgPuDF815nqfjhyKwo76014+8y8rbRRjMSrE1mZBxG9eLn7bX23BsR/6M+9kvA70TEzRHxA+A04LgJjfm3EfF3Uc2d/Rrwgjr9vwDHA9QvpY+r08yIKUWb/E+wO3fGjvM1fkT15iaHZwCflvRwI+0hqjkb353I+zTgmsb+t9s2dOK8uwC3PToQxuOAWxp5/mZtIyJ+VOd7AvDkjmXvioi7E206QdJbG2mbqK7bNBizoBgzg1vWz/IIbbTRjARrkxkZT6d6kTxtTu3T2PE56ttUz8T7NNL+prHdfLa7CPh/JT2NygEVVO4sM2JK0Sd33BbHtG/MLcC/iIgvZpS/jUqE1jigQ1tuoRrN2ishaEOWfbKk3SPi+1OO/VZE/FZGPT+sf/8Y8L/r7ae2bMvKUor4mGJJfUGb6T+k+vsFHrE0PaVxvI02mpFgbTIj4xbgAE0PiHQr1QuiNQ6gmgLyPWBm9O+I+L6kz1E5oZ4D/EGkv/w7aB0b6FllbJSiT7ZKLo7vAXuuTX6t+Y/Ab0l6BoCkp0g6OlH+U8CJkg6R9GPAv5m3IRFxG/A54EOSnlhP0v37kl62gLJ/DPwHSXtI2kXSz9aHfw/4VUkvUsXjJf0TST8+pZ47qN66/1+SdlIVqOXvZ9+AFWYeK1KOWEk6StJfS7pJ0qlTju8q6cL6+FckHVinHyHpGkl/Uf9+RaPMsZK+Luk6SR9opD9D0hX1sask7T9xrieqCnozOd/ArAbfo5ozMosbgN3qv/FdgHcDzfmsbbTRjIChtMmYDvxPqpfeZ9XPFLtJeml97A+AX5f0TElPAP5v4MIWL6T/C/BGqrlus2ySXwX+sarAbE8F/vVcV2I6UZI2ueO2ICLir6iE4mZVUROfBpwNXAx8TtLfAn9OFdhjWvk/Bv49cCVwU/27C2+ksiJeTzXx/yKquWhDl30D1bzAvwJupxaxiNhKNc/tI3WdN1FNHE7xy8DbgTupghN8KfP8K0/fAlSPdnyUav7hIcDxkg6ZyHYScHdEHEQ1p/H9dfp2qkAUPwmcQBV8BlWLIH8QODwinks1x/HwusxvAx+PiOcDZwDvmzjXvwP+DLOqvA94d61zb5uWISLuAd4CnEP1EuaHQDPSWrY2mvFQ0sORWX2immf/KqpAat+h0phj68PnUv2/+gLwv4B7gbdOqSbFxVQ2ye/F7PUjL6CaH/ctqpfeF7Y4h+mRUrRJY26cMWZHnv/858cf/dEftS73jGc845qI2DztmKSXUAWO+Ll6/zSAiHhfI89ldZ4v15O3/wZ4StMeUk/S3k41d+D5wPsiYi28+xuAl0TEWyRdB/xcRGyry9wTEU+s8x1K1SH/E2BzRJzS+mKNMQtnCG0yxpg+mEefxqpNHnEzZsUY4K32fuwYOGJbnTY1T20luYdq2YYmxwDXRhUF8Cbg2arWw9mZKrT72hzNr/FoKOVfAH5c0p6SHgd8iKrjZoxZMTziZowZK6Vok4OTGLNizCkoe0na2tjfEhFb6m1NyT95kpl5JD2Xyj55ZN3GuyW9mcoW8jCVlXVt3tPbqBa1P5HKpvJdqknhbwEujYhbGhFLjTErwpgfdowxG5tS9KlTx61eFPpsYCfgnIg4q5dWGWP6ZvuMIf9t7BixdH+qiFvT8myrR9CeBNwFUAcX+TTwxoj45lqBiLgEuKTOczJVOHeiWsD0F+v0JwDHRMQ9tWXzZyS9hSrs8iZJP4iIxwRLycH6ZIwZI9YmY8y8zN1xawQ0OILqoe5qSRdHxPV9Nc4YsyMDDeFfDRws6ZlUo1/HAa+fyHMxVfCRLwOvAa6MiJC0O/BZ4LSYCN0uae+IuF3SHlSjaa+r0/eiWs/vYapFT8+tr+2XGmVPpJrjNm+nzfpkzAIZu71oLFibjFk8JelTlxG3w4CbIuJmAEmfBI6mijQ4FUll3DVj+mN7RDxl/WyP0rf4RMSDkk4BLqN6A3xuRFwn6Qxga0RcDHwMuEDSTVQjbcfVxU+hith1uqTT67QjI+J24GxJL6jTzoiIG+rtlwPvq/XgC8Cv9XpBFa30aa+99ooDDzxwgGYYs5pcc801S9emQmn97GR9MmZHNrI+dem4TQto4HDNxrTj220LDCE+EXEpcOlE2nsa2/cCr51S7kzgzESdxyfSL6JaQmJWe84Dzlun2bNopU8HHnggW7duTR02ZsMhaRTaVCCtn52sT8bsyEbWpy4dt5yABmtzW07ucB5jTINSxGdg1tWnpjYdcMABi2iTMUVjbcqi9bOT9cmY7pSiT106bjkBDagj122BHa2Szahxze3HPe5x627vvPOOzU7lS9XbJNWOVFS7nGh3uRHxml+i5vZDDz20bp5UPV3S583XR9l5zjVE+3LTU59Lc/vhhx9eN888lCI+A7OuPjW1afPmzb3c1NKjYfq7Z2bh70cWrZ+dUvrUvN99PpvklE+dO/UdyMnTtUzq/2zOdbfNP1RdfbajCzmf46Lb0efnssp06bjlBDQwxvRISRNsB8b6ZMwCsTZlY20yZsGUpE9zd9xSAQ16a5kxZiqliM+QWJ+MWTzWpvWxNhmzHErRp07ruE0LaJBLc8hzp512emQ7ZYlsbm/atGmHunbZZZepdaXqTdkm2+bpYqeE9LD+gw8++Mh2ym7XtFO2rX/Wl7etpW8R9s0uls+csqk8qXs/eay53fxcUtu2Si6GLvpkjGmPtSmPLtqUY1HMfQZpa31s+794HntkTr05tL22Pi2Aqbpyzte2fZP3qMu5F23xXPTnUoo+deq4GWMWTyniY4wpC2uTMWaslKJP7rgZs2KUIj7GmLKwNhljxkop+rS0jlsq4mOOPXK33Xbboa6mVTJlm0ydIyeKZRdrZS4pq2TKEpmy6rXNM3nuVFTEHNraD1NlZx3LaVNba0fO9U+et/m5ND+vHAtC2/s6WU8p4mOMKQdr02LoYqWbzNNlGkKXCJPzkGsVbNOmISyDs2jbjj4jheacu8u1zmPl7DrdqA0l6ZNH3IxZMUoRH2NMWVibjDFjpRR9csfNmBWjFPExxpSFtckYM1ZK0adRWCXb2iabdkjY0UaZslem6s2JQpkTATMn2uRkviYpW17Ogs456bMiJKbsmH1ZH3NsgpP1ty3TxU7ZJBX9cfIeNe2RDzzwwCPbKTtCThTQXEoRH2NMWVibhqcvC2Bu+T4jRHehbb1dojZ2jYrYNn9f6bnnW6RFdEwRJkvRJ4+4GbNilCI+xpiysDYZY8ZKKfrULaKGMWahrE2wbftjjDFDsgxtkvRkSZdLurH+vceUPD8l6cuSrpP0dUnHdjqpMWblWMZz01D6NAqrZE40x5RFcfJYKmJkyjbZdvHueaySORbMJm0jIebY+2ZFSEzZKNtGXsxpa47FMzdfl3Y0SZVt2iGb25P7zc80df+7Rh01xhjzGE4FroiIsySdWu+/cyLPj4A3RsSNkp4GXCPpsoj4/qIb2wddIhPOiirZl3Uvde4U89j+2jJEFMVZ5XPO18UCOI89sm37+mKWPbdL+oowiD7ZKmnMiuERNGPMGFmCNh0NvLzePh+4iokHo4i4obF9q6TbgacAK9lxM8bMRyn65I6bMSuGO27GmDEypzbtJWlrY39LRGzJLLtPRNxWn/s2SXvPyizpMGAT8M15GmqMWV2W8Ow0iD6NouPW1jaZS6pMjt0xJwpls92pPJP5utgxUzStfqmFoZvbuVbJnEiIXSyRs6ySKWtnzoLkbSNMpphllbz//vunlknd/1T+eXDHbbyU/tn0ZVsp/T71RZ82oUXc8znPsT0iNqcOSvo88NQph97V5iSS9gUuAE6IiPn/MYyItt+Pyf+fbcv3Ffkv13LZtn2paQs5NM/V5blhsq6+8uTaUftctLuPerou+r7kqJLrvlRahj6NouNmjMnHD73GmDEyhDZFxCtTxyR9T9K+9dvsfYHbE/meCHwWeHdE/HnvjTTGjJ459GnmS6W6zoXrkyMmGLNCLCNymzHGrMeStOli4IR6+wTgM5MZJG0CPg18PCL+a9cTGmNWjyU9Nw2iT+64GbNiuONmjBkjS9Cms4AjJN0IHFHvI2mzpHPqPK8DfhY4UdJX65+f6npiY8xqsYTnpkH0aRRWyVQ40q5hSlNz51Lpqflnqe3U0gPN9Nzyze3UEgU5oV2bHvb77rvvke1Z88RS88Zy5q/1Ncdt0lOes5RBajunTU1S3vbmHLUHHnhgatnJMs18XedopnBHzBgzRhatTRFxJ3D4lPStwJvq7U8An1how5ZAzlyg3OUAUrSd15bKM1To9y7fv3nKDjG3r9Cw+FNZ5Hw8KEefRtFxM8bk446bMWaMWJuMMWOlFH1yx82YFaMU8THGlIW1yRgzVkrRp6V13HKGg1Pbk7azHEtkTv4c22TK3piyTULa+thM37Rp09T05nbKbpeySu66666PbM+zHECOzTBlfUzVn7IiTtafskq2Tc+5nhyr5OS9b5Zv5kt9D/oMo16K+BhjysHatFy6TifpYpvsq02T35+c9rX9znUpO6uuLufuYrPsWqYvy+s8n0/O/Z/nulP1lKJPHnEzZsUoRXyMMWVhbTLGjJVS9MlRJY1ZMYaI3CbpKEl/LekmSadOOb6rpAvr41+RdGCdfoSkayT9Rf37FY0yx0r6uqTrJH2gkf4MSVfUx66StH+d/lOSvlzn/7qkY3u4XcaYBeGIt8aYsVKKNo1ixK1tFJ3J9JTdMcf6mMqTY4nMiQqZm69plWxaHFNWydRQdCrC4SyrZI7dMYeciJE5VsfJ/bbbqXpzrJ/N/M37PZm/eT9TkT+HsEpOa0tXJO0EfJQqXO024GpJF0fE9Y1sJwF3R8RBko4D3g8cC2wHXhURt0p6HnAZsJ+kPYEPAodGxB2Szpd0eERcAfw21Zol59cdvfcBbwB+BLwxIm6U9DTgGkmXRcT3e71gY8wgjPlhx8ymi5UuRVf7ZV82uVSdQ5ETwXooe2TqfG0jfOZYH3Om0OS2KSd/V0rRp3VH3CSdK+l2SX/ZSHuypMsl3Vj/3mPYZhpj1hjgrfZhwE0RcXNE3A98Ejh6Is/RwPn19kXA4ZIUEddGxK11+nXAbpJ2BZ4F3BARd9THPg8cU28fAlxRb//p2rki4oaIuLHevhW4HXjKrIZbn4wZDx5xexRrkzHjohRtyrFKngccNZF2KnBFRBxM9QD2GGuVMaZ/5nkwyhCg/YBbGvvb6rSpeSLiQeAeYM+JPMcA10bEfcBNwLMlHShpZ+DVwNPrfF/j0U7cLwA/Xo/QPYKkw4BNwDfXaft5WJ+MWToDadMqcx7WJmNGQUnatK5VMiK+sDafpcHRwMvr7fOBq4B39tGgnEiSk0OnOVEi+1poO2V1zLVKpiJGNq2SqQiTKYtnyqLYPG/KMjhtf40uw/opu2Jqu2k9nHWsmd60grY9XyryZPO+NvPMsucuwh7ZZE5B2UvS1sb+lojYUm9Pa+jkSWbmkfRcKvvkkXUb75b0ZuBC4GHgS1SjcABvAz4i6UTgC8B3gQcbde0LXACcEBEzvbuL1idjTJoxP+wsmmVq0zwWu7ZWuhzaPjfMKj9E+4airzb1eW1dImC2zT/r+5ejEUNYZHPPvQrMO8dtn4i4DSAibpO0d49tMsb0z/aI2Jw4to1HR8MA9gduTeTZVo+gPQm4C6AOLvJpqvlpj4yQRcQlwCV1npOBh+r0W4FfrNOfABwTEffU+08EPgu8OyL+fM5rtT4ZY8aItckY04nBo0pKOlnS1om3/caYORlgyP9q4GBJz5S0CTgOuHgiz8XACfX2a4ArIyIk7U7V0TotIr7YLLD2UFLP43gLcE69v5ekNe05DTi3Tt9E1QH8eET819Y3piVNbbrjjjvWL2CMmUlJdqRlY30ypl9K0aZ5R9y+J2nf+o3RvlRBBKZS27G2AEhqWquYd3vWAtypfG0jSabSc+yRzfTJY20X3c6xSjZpftma7ciNFpm6l6nIialzp2yJTdtjc/v+++/foa6UVbLLdqrOVATMlNV2FouwbfQtKBHxoKRTqCJC7gScGxHXSToD2BoRFwMfAy6QdBPVSNtxdfFTgIOA0yWdXqcdGRG3A2dLekGddkZE3FBvvxx4X60HXwB+rU5/HfCzwJ61jRLgxIj4astLytKnpjZt3rx5vCptzIow5oedkTDXs1Nbfcr5Hz2Zb1mWw9zzDt2+rvWn7H1tbX9dbYI55bvkadumWd+5nGk6fVKKPs3bcVt7+35W/fszvbXIGDOTIcQnIi4FLp1Ie09j+17gtVPKnQmcmajz+ET6RVSRKSfTPwF8olXDp2N9MmYJlPJgNCDWJmOWRCn6tG7HTdIfUL0h30vSNuDfUInOpySdBHyHKQ90xpj+GfsQ/qKxPhkzDqxNO2JtMmY8lKRPOVElp741Bw7vqxFdrZJ9La7dzJOyPubYI5sWyNwyqXbkRMxMDSs365w1NJ5jL03dv9Swd8oe2bRENrcno0o2I0amFhLvkqfZ7mZ6ikmraSoq5SKEoRTx6YNF6JN5lDF+9/qy1Yzx2sbYplmsWnuHZNHaNI9Vr61Vbahof32R+l/c1/XMsv2lyJ1mMa3OeZ4tcto01P2Yt55FUYo+zWuVNMYsiVLExxhTFtYmY8xYKUWf3HEzZsUoRXyMMWVhbTLGjJVS9Gl0Hbe2i2nPKpNjm2wbVTK13bRH5kaVzIkemWpTjm0yZ3tyP+f6cqItNq2ETStiyio5aVdMWRxTdaXScyylzbY27/EslikApYiPMaYsrE3LI/X/rc+ojV0iCg4VITHFEIt3z7KajtGyPYRdNFW2yVCLZvcZ+XOVGV3HzRiTpqQJtsaYcrA2GWPGSkn65I6bMStGKeJjjCkLa5MxZqyUok9L67h1iSSZG1Uyx3LYZaHtnDyzyuREmMyxSqbsfTnWysl8qciaKStojlWyGeXx3nvvnVrn5ALcOfbIZvua6anrTlknmgtz5yxyPllXDkNZB4wxZixYm4any0LP80RF7MuumEPuYs1DLxqdS9t25FgL+/wbamsXbftZDxV9dKiolKXok0fcjFkxShEfY0xZWJuMMWOlFH1yx82YFaMU8THGlIW1yRgzVkrRp1F03NraI2dZ/dpGj0zZFVOWwdR20/Y3aZVsa8HMsXXmXHPqnk1aK5vHcq4pFQGzScoq2by2lO1xskyOPTJl32xr5cy9ZzkRvIawZ5Q0wdYYUw7WpsXQNnpkV9vZUP/HptU/61xDRL0cirbWx7bWyq5tyjl3Dm3tl7nn8LPTbEbRcTPG5FOK+BhjysLaZIwZK6XokztuxqwYpYiPMaYsrE3GmLFSij6NLqpkKgrgLPtbc7/tItqpPDk2wZztyf3m+VKRGlO2vxwLZcrel1pke9a5m+1LtTV1jlTUxpxzwY6RJFPRJ++7776pdaW+K802pRYIb35WuYu+L9qeUYr4GGPKwto0PuaxrXWJXNm2TYugr3bPY/tLkbqvi/4byvl821pb57nfQ0WSTJ1jlUnHOzfGGGOMMcYYMwpslTRmxSjlrZExpiysTcaYsVKKPo2i49Z2Ae7JCH9t7ZE5C2rnRH9M2SknbX850SNzrJw52ylLZCr/rHPn3JuUtbVJ848lZY9sRnac3G/aI1ORJHPa1GxHs/5mO5oLhM9aaDxlo3RUSWPMRsTatHi6Lpjc1gLXl81w0VbJvs43q562ERbH+LfS5XMfY/TSJiXp0yg6bsaYfEoRH2NMWVibjDFjpRR9csfNmBWjFPExxpSFtckYM1ZK0afRRZVsbudES5zc72KPbGuJzF2AO2fx6lR6TgTHLtbKyfI5lppwUoAAACAASURBVM2cyJWpYe/mdaYiT8KOkR6b202rZMrimGprk5wFvlPXNlmvo0oaszz6+nsYy+K9TVbtb33V2rvqdLWwDWGBm2dx7S7n7itPl/yz8jUjWHdhnnak/h77yjN09NG+KUWfPOJmzIpRivgYY8rC2mSMGSul6JM7bsasECVNsDXGlIO1yRgzVkrSp1F03HIW4J4V4a9LVMS2USLbbk/W1cWW2Na6OI9Vsq3tNBVdsW2kxUk7QfOeNa2SOdbRVL1Ne+Suu+469Vyp+zS5aHnbBeT7tAuUIj7GmLKwNg1Pnza0Lla3Pi2RKeaJlDlvnlT+RSwmvYjPNJWnSdvFuNuet48yXShFn0bRcTPG5FOK+BhjysLaZIwZK6Xo0/TIDQ0kPV3Sn0r6hqTrJP2rOv3Jki6XdGP9e4/hm2uMWRvyb/NTItYmY8aFtelRrE/GjItStGndjhvwIPCbEfEc4MXAr0k6BDgVuCIiDgauqPeNMQMzxMORpKMk/bWkmyQ95m9Z0q6SLqyPf0XSgXX6EZKukfQX9e9XNMocK+nr9UPLBxrpz5B0RX3sKkn7N46dUD/Q3CjphHWabW0yZkS447YD1idjRkQp2rSuVTIibgNuq7f/VtI3gP2Ao4GX19nOB64C3pl74pzlAFLzmSbD7XeZ19Yl7H9zntSs5QCax3LmirWdy5aaW5azTMCsdqTm5uWcO7U9ee41Jue45cyXa/5hNcs3lxZIhfrvMtcwt31DeLaHEBRJOwEfBY4AtgFXS7o4Iq5vZDsJuDsiDpJ0HPB+4FhgO/CqiLhV0vOAy4D9JO0JfBA4NCLukHS+pMMj4grgt4GPR8T5dUfvfcAbJD0Z+DfAZiCAa+p23J24F4NokzGmPWN/2Fk0i9aneeacdZnLtcww78ucQ5ZTb5c5dTl5cv/O2v49dpm/Ns99XeR3qCR9yhlxe4T6LfsLga8A+9TCtCZQe/fdOGPMYxngrfZhwE0RcXNE3A98kurhosnRVA8ZABcBh0tSRFwbEbfW6dcBu0naFXgWcENE3FEf+zxwTL19CNWbZoA/bZzr54DLI+KuurN2OXBUzj2xNhmzfDziNh3rkzHLpxRtyu64SXoC8P8B/zoi/neLcidL2ipp6zwNNMbsyAAPR/sBtzT2t9VpU/NExIPAPcCeE3mOAa6NiPuAm4BnSzpQ0s7Aq4Gn1/m+xqOduF8Afrweoctpx2PoQ5vuuOOO9QsYY2bijttjsT4ZMw5K0aasqJKSdqESnv8cEX9YJ39P0r4RcZukfYHbp5WNiC3AlrqeqXeiOVyaYw2cDM2essC1DenfZQmAlG1yst6U5TNlZcyx5LVdAmDSrpizBEBOXW2tkimr4yya+XLal5OnrQUV8qy+QzGnoOw18fJkS/23CTCt0ZMnmZlH0nOp7JNH1m28W9KbgQuBh4EvUY3CAbwN+IikE4EvAN+lmg+S044dG9WTNm3evHm8Km3MijDmh51lMIQ+tf0fk2sTXKb1MYec71aOdTHnfsxzz4ayY06rfxY5bWqbPk87UvmH+v7O05ZVJSeqpICPAd+IiN9pHLoYWAsecALwmf6bZ4xpMs8b7VqstkfE5sbPlka123h0NAxgf+BWduSRPPUI2pOAu+r9/YFPA2+MiG822npJRLwoIl4C/DVwY51+a0T8YkS8EHhXnXZPZjsewdpkzHjooE1FYn0yZjyUpE05VsmXAm8AXiHpq/XPPwbOAo6QdCNVUIOzBmynMWY4rgYOlvRMSZuA46geLpo0HzZeA1wZESFpd+CzwGkR8cVmAUl717/3AN4CnFPv7yVpTXtOA86tty8DjpS0R13myDothbXJGDNWrE/GmN7JiSr5P5huYQI4fN4T51jNUhESZ1klc2yQfVkic7ZnnXuWfXFaesoqmWP1y40qmbIW5kTDzLF1NtNnvdFoG7Wxy/Y8USFTZZoMZaHs+01QRDwo6RSqTtJOwLkRcZ2kM4CtEXEx1ZvjCyTdRDXSdlxd/BTgIOB0SafXaUdGxO3A2ZJeUKedERE31NsvB95XW6e/APxa3Y67JP07qo7kWpm7ZrR7EG0yxszHot9S15FoLwQOBL4FvC4SUWglPRH4BvDpiDhl6LYtQp9ybGSzohG2tc+lyrZlHvtbl8iLi5jOMBaraVcr43rpba2VXe/LmJ+d1mMofcqa42aMGQ9DiE9EXApcOpH2nsb2vcBrp5Q7EzgzUefxifSLqCJTTjt2Lo+OwBljVogl2IvW1kQ7S9X6k6eSDq3/74A/W1jLjDGjohR9arUcgDFm+ZTk1TbGlMMStKm5TMn5VNFrH4OkQ4F9gM91PaExZjVZwnPTIPq0tBG3HKtfyrY3GbUxx8q42267TU1va5VsWh1zzjWr/KwFntdoa/VLWfhS93LyWGo7FblzVuTF9Zj1h9E81owk2dxu5kltp/KnyF04PHWtbe/BPLgjZowZI0vQph3WRFubV9uknk/7Iar5ZkVZqNvayOaJKjmEBbBrxMK2tsm2DGV7zLGjLuLcTdraYttewzz23KEoRZ9slTRmhfAImjFmjHTQpllLlSDp88BTp5R7V2b9bwEujYhbxjIPyRizWObUp5naBMvRJ3fcjFkx3HEzxoyRObVpe0RsnlHnK1PHJOWsifYS4GckvQV4ArBJ0g8i4tR5GmuMWU3m0KeZ2lTXuXB9GkXHLWUBTNkmZ1klmzbF1HbT1piK+JiKPJmyRz7+8Y+feq7J8rMsi2vkfLlSPfM+o0qm7kEzPSe6YtOu+OCDD66bPrmf2n7ggQfW3W7mf+ihh6aeO8WsSFSp7+kiole54zZexhgBa4xt6pO+rm+M17ZqLOEeri1TchaJNdEi4pfWtiWdCGwupdPW1mqWa5XMqTfnf2gOfUYdbNumtvdvlu0vRWr6RF/RH3PLDPE80tX6mWOZXfH/X4Pok4OTGLNiODiJMWaMLEGbpq6JJmmzpHO6Vm6MKYclPDcNok+jGHEzxuTjjpgxZowsWpsi4k6mTOiPiK3Am6aknwecN3jDjDGjoxR9GkVUyZzFkGctwN0lMmTb7baLcc9qR2oIPRVFMUXOvZxllUwdy4kq2cXuOes6m7bGlD0yZaFslk1tp+5rTrTOyf22C353wSNoxpgxYm1aPIuI0Df0Z5oT1XAWqTI596PPxaRzoyquR1vL5mT+RUau7DMS6SJsnaXok0fcjFkxShEfY0xZWJuMMWOlFH1yx82YFaMU8THGlIW1yRgzVkrRp9F13FLR+mYthpyy7rVNz7EMNtNz7JSQjmKZskqmLH1to02mrH6zbH+pe56Tp0nKiti8huZ285on91P3I5Xe5Z7l2E5n5VsEpYiPMaYsrE3Dk7KUzbOo9TwLKLdpUw5dbX9daGsrnBVVsq0ds22erpEaF/mcMs/3odAFuAdhdB03Y8xsShEfY0xZWJuMMWOlFH1yx82YFaKkCbbGmHKwNhljxkpJ+jSKjlvb4eZZiyHnWN1yttvaJlMRGCf3U+XbDqHn2CW6RkjMWVw7h1S7Z0WVTB1ru932nuXaS3Mstn3dv0lKER9jTFlYm4an7f+S3EiIOel9tSnFLCtizjNS2za1LbvoaImLnobRhUVENe1abyn65AW4jTHGGGOMMWbkjGLEzRiTTylvjYwxZWFtMsaMlVL0aWU6binL27T9nDJrtF08OcdOOWuB66ZVsllXykLYPHfKApgaPs5ZJHq9Y9PIua9d7I2T+22jR3axR6Y+q8lF31Of9yKiTZYiPsaYsrA2jZuuizXn/G9dhLVtkREmZ7WnSwTNFG0/h7HYKcdqj0zVu8qsTMfNGFNRivgYY8rC2mSMGSul6JM7bsasECVFRjLGlIO1yRgzVkrSp5XpuKWsc5P7Oba6oRdonrVYc459sWm9a15PylqZ0+6c9Ml6Z93zaXmapD6TnM9qcj8nEmXbiFNt7ZGTVsmU5bWZPpRtshTxMcaUhbVpeLrY52Z9Pm0jTC5ioee2kSH7iuCYu8j5EJa+nDqHim7ZV56u7fA0k9msG1VS0m6S/qekr0m6TtK/rdOfKekrkm6UdKGkTcM31xiz9uaozU+JWJuMGRfWpkexPhkzLkrRppzlAO4DXhERLwB+CjhK0ouB9wMfjoiDgbuBk4ZrpjFmDT8cPYK1yZgRYW3aAeuTMSOiFG1a1yoZVet/UO/uUv8E8Arg9XX6+cB7gd/t2qBU1KLm9qRt74EHHnhk+/77739kO7XYdWqR5E2bpr/4GuoDHCJKUNuokJP7KVvigw8++Mh2yr6ZsjE28zfryY0qmWpf6hpSpOyRzQXSm9+BXXfddWr65H5zO2WZtVWyfxatTZltWsRpTIOS7/mq6UbJn0VbhtKnLvbIVft82ra3r7+X3Hrmed7qm8n6c+yHuVbQ9eqxPi2HrAW4Je0k6avA7cDlwDeB70fE2lP4NmC/YZpojFljnjfapYjVNKxNxowDa9NjsT4ZMw5K0qas4CQR8RDwU5J2Bz4NPGdatmllJZ0MnDx3C40xOzBmQVk0fWnTAQccMFgbjdkoWJt2xPpkzHgoRZ9aRZWMiO9Lugp4MbC7pJ3rN0f7A7cmymwBtgBImnrXUjczx3o3uZ+y9KW2UxbMVJ6225P0Zfvrk5Q9MmfIPWVp7DOqZNsooKm2pmyMqaiSqfTJ/VS9k9FF+2KI74eko4CzgZ2AcyLirInjuwIfBw4F7gSOjYhvSToCOAvYBNwPvD0irqzLHAu8q67zsxHxjjr9ACqL0O71sVMj4lJJuwDnAD9NpU0fj4j35bS/qzZt3ry5DEU3ZomU8mDUN4vWp3k+h76mawzFEItupxjqXrStNyd/bpTwHIaYsjPLypnzPLxqdsxFkBNV8in12yIk/T3glcA3gD8FXlNnOwH4zFCNNMY8St9D/pJ2Aj4K/DxwCHC8pEMmsp0E3B0RBwEfpppgD7AdeFVE/CSVDlxQ17kn8EHg8Ih4LrCPpMPrMu8GPhURLwSOA/5Dnf5aYNe6rkOBX5F04Ix2W5uMGREl2ZG6Yn0yZlyUok05I277AufXD3ePo3rg+iNJ1wOflHQmcC3wsQHbaYypGUBQDgNuioibASR9EjgauL6R52iqSfQAFwEfkaSIuLaR5zpgt3p07lnADRFxR33s88AxwBVU1qAn1ulP4tE3zgE8XtLOwN+jGsH73zPabW0yZkSM+WFnCVifjBkRpehTTlTJrwMvnJJ+M9UDX6/k2AcnbXUp+10z2mTTQpmyU7Zd9Dlne71j0+hqx5zGrKHrVL2ptuZEr+qyKPrkfs61piIdpbZz7I2pSKSz8qXskX0N9w/0Jmg/4JbG/jbgRak8EfGgpHuAPalG3NY4Brg2Iu6TdBPw7HrEbBvwaio7JVQdwM9JeivweKo30VB1CI8GbgN+DPj1iLgr1ehFa5MxJs3Y31IvmkXoU9fI1Iu0R85j++vyzNPl2oa6F20/r67t6MuO2ZZ56lzE968UfWo1x80Ys3zmFJ+9JG1t7G+p51AATFPMyZPMzCPpuVT2ySPrNt4t6c3AhcDDwJeoRuEAjgfOi4gPSXoJcIGk51E9zDwEPA3YA/jvkj6/NhJojBk3pTwYGWPKoxR9csfNmBVjTvHZHhGbE8e2AU9v7E+bML+WZ1ttZXwScBeApP2pIqa9MSK+2WjnJcAldZ6TqTplUM2XO6rO82VJuwF7Ua1t9CcR8QBwu6QvApsBd9yMWQFKeTAyxpRHKfo0TNg7Y8wqcTVwsKRnStpEFTDk4ok8F1NNpIdqYv2VERH15PvPAqdFxBebBSTtXf/eA3gLVcRIgO8Ah9fHngPsBtxRp79CFY+nisD2V71eqTHGGGPMijKKEbecldhneZ1TSwWk0rtsp+bKpbZhxzlQzeubnDe13jX09bZgcj5ZTnjWnGUC+pwX2GVuX868trZz4sZE32+N6jlrpwCXUYXnPzcirpN0BrA1Ii6mmkB/QT137S6qzh3AKcBBwOmSTq/TjoyI24GzJb2gTjsjIm6ot38T+D1Jv05ltzyx7gR+FPh94C+prJm/X88TMcasAKW80R4zqflDbecz5Zbpi5x5XLO+P22XA1jm/+/UZ9TXfLI+5zO2bUeXuYazzjf0PZvWllVlFB03Y0w+Q4hPRFwKXDqR9p7G9r1U4fony50JnJmo8/hE+vXAS6ek/2DaOYwxq0EpD0bGmPIoRZ/ccTNmhSgpMpIxphysTcaYsVKSPi2t49bFnjbLKpnazrE+5lgiU3maSw80tyEdIj5lP2y7fECTHJvgJM1rSjGr/Bo5NstcC2SX5Q5y7JEp5mlfzhIKfVKK+BhjysLaNDxd7GJjtA92LT+QAyXrXG3tn32F559VT1tLZOr+tT137jX3db55KEWfPOJmzIpRivgYY8rC2mSMGSul6JM7bsasGKWIjzGmLKxNxpixUoo+jcIq2bSzpbZzh4abVrWUxbGLbTJlibz//vsf2U5Fi5xs3847P3r7+4qYlDM0Pml77GI7aGsJSNlAc62IKZrXlIri2TYaZup7Mrmfske2vYZcShEfY0xZWJs2Jl3sffNEksz5nnWxZvYZ1bBtFMWuUR7nmW4077lTdebaSxetF6Xok0fcjFkhSppga4wpB2uTMWaslKRP7rgZs2KUIj7GmLKwNhljxkop+rS0jlvTztbcblremlbCWVENcyIYtl1ou2mD3GWXXaamz2OPbNouU9fXZaHoPheQTrUj9Vk001N/ILnRGHMWHk/d85QdtXnvm21NfT6zIoWmoo4OsWD6JKWIz1jpyy4yFsbYpj7pK+pYn/dpjG1aBKvW3lWki+1vqKjHbe2AqTy55dtaM/uyTc5DX1FAU88W80SYbD479TU9ZhVYtfam8IibMStGKeJjjCkLa5MxZqyUok/uuBmzYpQiPsaYsrA2GWPGSin6NIqokqkh3FzbX8p+l4r8l7PQdsoy17TeNa2S81g5m+1oa5tsm55qzyxyFs5ue47cqI250SfXIye6ZY5dtvl9mNxPWSWHsE2WNMHWGFMO1qblkhOBcZ4IhItkngWuc2hr5ewaSbItbT+73DYNYWtc1b/xkvTJI27GrBiliI8xpiysTcaYsVKKPrnjZsyKUYr4GGPKwtpkjBkrpejTSnbcZkUgzLFH5izAnWOV7LpYeKqtzXOkIm7m1J+K8phr68yJJJkq2zaK56wFrmflWyNlS8yxLqbqnyfq5VBRu5qUIj7GmLKwNi2PLotgD0VXS2PbaRlDR57MtXK2XeA65++mz+iWGymSZJNVbnuT9Z/GjTHGGGOMMcYslZUccTNmI1PKWyNjTFlYm4wxY6UUfRpdxy3HSjhpR2trj8yJJJljm2zaB1PRMCfpa3HflPUxdZ9SeSaPpRa1Ti0E2UxPfQ6phaxnLXDd3M+Jzpj6fHMslzlW27FYJUuKjGSMKQdr03gY6nPoKwpj16iNXWyhqbI57ZgnmmPq3F3qmVVn2+vI+UwX/Xfd1mqaW2cp+pRtlZS0k6RrJf1Rvf9MSV+RdKOkCyVtGq6Zxpg11gSozU/JWJuMGQfWpsdifTJmHJSiTW3muP0r4BuN/fcDH46Ig4G7gZP6bJgxZjp+OHoM1iZjRsCitUnSkyVdXneCLpe0RyLfAZI+J+kbkq6XdGCnE7fD+mTMCFj0c9NQ+pRllZS0P/BPgN8CfkPVeOUrgNfXWc4H3gv8btbVkB6GzYn8N2nna2uPbEZtbGuPTEV5nBVVsu3QdU6EyhzrYpNZX8JUXanPIhVhMvU5NO9fc9Hy5vakVTJnUevUuVMWzJSFMsf2OOv+pY4NZTXYAB2xbIbQJmPMfCxBm04FroiIsySdWu+/c0q+jwO/FRGXS3oCMHz4X4bRp3nsczn5ulgF+1oQu2u9TXIskV2tn10snzmRGlPPEG0tnrPqzUlfNENFPy1Fn3JH3P498I5GZXsC34+ItSfhbcB+mXUZYzrgEbcdsDYZMxKWoE1HU3V+qH+/ejKDpEOAnSPi8rqNP4iIH3U9cSbWJ2NGwhKemwbRp3U7bpL+KXB7RFzTTJ6SdepVSjpZ0lZJW9c7lzFmNvM8GJXacetTm+64445B2mjMRmFJ2rRPRNxWn/82YO8peX4C+L6kP6znmn1Q0vQoXD1ifTJmPCzpuWkQfcqxSr4U+GeS/jGwG/BEqrdIu0vauX5ztD9w67TCEbEF2AIg6ZE7kbLV5VgXJ616TftiytaYskem0lP1NNNT9sbJ9qWONc/XPEdzu+2Q9jxD4KnIkKm2poaxc6JK3nfffVPTJ62SOREgU7bGZlubdszUuVP15zLUsH6KUjtic9CbNm3evNk31ZiOzKlNe0282N1S/20CIOnzwFOnlHtXZv07Az8DvBD4DnAhcCLwsXka24LB9Sl1v3NseJP5VomcKSE51seu1s++7JE5n9eqflZd6WJHnVVXJjO1qW7TwvVp3Y5bRJwGnFY38OXA2yLilyT9V+A1wCeBE4DPZDbSGNMBd9wqrE3GjIs5tWl7RGyeUecrU8ckfU/SvhFxm6R9gdunZNsGXBsRN9dl/hvwYgbuuFmfjBkXc+jTTG2q61y4PrWJKjnJO6km295E5dse+u2VMQbPccvA2mTMEliCNl1M1fmBdCfoamAPSU+p918BXN/1xB2wPhmzBJbw3DSIPrVagDsirgKuqrdvBg5rU75JTtTBpo0xZUuEtM2wbVTJVMTInEiSqYW5Z7W9ud1sU5cvTMo+mBshsUmOpSB1vhz7a1erZE4E0ZyIls2yTWZ953Iifw7FEB0xSUcBZwM7AedExFkTx3elinx0KHAncGxEfEvSEcBZwCbgfuDtEXFlXeZYKrvATsBnI+IddfoBVBN1d6+PnRoRl9bHng/8Jypb0cPAP4yIe9drf5/atAE7utn09V3v8x6P8fPqq019assi7tMSPouzgE9JOonKZvRaAEmbgV+NiDdFxEOS3gZcUUd1vAb4vUU2cqhnp7aWwXXa2LrMMuqcRSridZO20RxnpbeN8Lno+5E6X9tI5Kk8zfvd9W9/ERbRUvSpVcfNGLNchhhBqyfCfhQ4gmrY/mpJF0dE863PScDdEXGQpOOo1iI6FtgOvCoibpX0POAyYD9JewIfBA6NiDsknS/p8Ii4Ang38KmI+N06otKlwIGSdgY+AbwhIr5W17Fjj94YM0qWMbofEXcCh09J3wq8qbF/OfD8BTbNGDMiStKnLlZJY8wSGMCOdBhwU0TcHBH3U829OHoiTzOs7UXA4ZIUEddGxNrk+uuA3erRuWcBN0TEWji0zwPHrF0C1YgawJN4dHL+kcDXI+Jr9XXeGRHTh0ONMaPDNm5jzFgpRZuWNuKWigLYdkFsSFslcxbaTpVtbqciWqbyz7LVpaJVphaHTkV5bJIa+s9ZcHqSnC9raui/rVUylWdWvpzt1HclFW2y7cLcy2aAtuwH3NLY3wa8KJUnIh6UdA/V/IztjTzHUE2yva+ev/FsSQfW9b2ayk4J1YKzn5P0VuDxwNrk3p8AQtJlwFOAT0bEB/q4QGPM8IxJJ0slxzqWE71wvWPrkVM2tx1daPuda2vJy2132/uRQ9d71qUdbW2gqbJj0oQxtaULtkoas2LMKT6zwtrmrC00M4+k51LZJ4+s23i3pDdThbZ9GPgS1SgcwPHAeRHxIUkvAS6obZY7A/8n8A+BH1F5vq+p7ZXGmJFTyoORMaY8StEnd9yMWTHmFJ9ZYW23AU9v7E9bW2gtz7Z6LtqTgLsAJO0PfBp4Y0R8s9HOS4BL6jwnA2vDmScBR9V5vixpN2Cv+hx/FhHb6zKXAj8NuONmzApQyoORMaY8StGnUVgluyzGPVkmVT5lS8yxO6bSU5EkJ4eYU5EkcyyeqXpT6SkLauoeQdr6mLOd+uzaRn+c5zNtW2/b7Vn3LMdGOYRIDOS9vho4WNIzge8CxwGvn8izFtb2y1RrEF0ZESFpd+CzwGkR8cVmAUl7R8TtkvYA3gK8rj70HaoJu+dJeg7V4rR3UAU2eYekH6OKUPky4MN9X6wxpn/GPi+kdFK2ulmWt7YLQqfq7WIT7Grl7MummXrGybUodrGFdrnHk39zfUYanbeeMepASfrkETdjNjj1nLVTqDpOOwHnRsR1ks4AtkbExVRrDV1Qz127i6pzB3AKcBBwuqTT67QjI+J24GxJL6jTzoiIG+rt3wR+T9KvU9ktT4xKUe+W9DtUHckALo2Izw557cYYY4wxq4I7bsasGAON5F1KFZa/mfaexva91GuQTOQ5EzgzUefxifTrgZcmjn2CakkAY8yKUcobbWNMeZSiT+64GbNilCI+xpiysDYZY8ZKKfq0tI5b2zDyqTlxsONcs+b8sJw5UKklA5pzzlJz2VLbuXPcUtvNNqW8zDlz3JrtnrUcQOrepj6jVHrOvLac+WqTx1Lfg5ylBXLmvjXrbC4TkFpKYNY5cu5lV0oRH2NMWVibFstQIfab9BU6Pqds7nyteepaLz2nnnnmu/WVfxFt6rLEQy7LXCqgFH3yiJsxK0Yp4mOMKQtrkzFmrJSiT+64GbNClBQZyRhTDtYmY8xYKUmfRtFxy7Hk5YZmT9kg2y4Z0NYeOcsqmQrvmhqKbtaVsh/2uRxAWxtkqk05n1fO55CbLye9+Tmm7k3KXnrfffdN3Z7Mlzp3KrxwV0oRH2NMWVibFktf9rd56mpbNqfOrnbPtt+/lG0vt02LsKpOI9de2rZ9XZZ4yF2uoO1ntMzvx1gZRcfNGJNPKeJjjCkLa5MxZqyUok/uuBmzYpQiPsaYsrA2GWPGSin6NLqokm3tdpP7OREqm5bIVHrTitjM0xyqbeaZhxwrXcr2l7JKNuvJjcqZY33MsUq2jSqZ+5n2FT0y59zNe3zvvfdO3YYdrZPNSJQ5tsmulCI+xpiysDaNjz4tfG3r6hqlsK+okm0tkV0jNeZcd1sL4Dxt7XLueWyQCT1b5gAAG6FJREFU66XPU77P728p+uQRN2NWiJIm2BpjysHaZIwZKyXpkztuxqwYpYiPMaYsrE3GmLFSij6NwirZZaFn6GbLaxvxsGlRbNrqmuRGIcpZhLx5jpSVMxWFssmsxaBzoke2XWg7J08q4uOsMjnWxxyrZCoqZCqq5N/93d/t0L7mfjNf7qLnXShFfMzq4e/eYlm1+71q7V1FcmxrQ0XfG2Jh6S75J0k953SJmDkPbS2KQ32mOfbKLucY6u990d/rVcMjbsasGKWIjzGmLKxNxpixUoo+ueNmzApRkk/bGFMO1iZjzFgpSZ9G0XHLsQ+mrHeTx1ILcOdY91KWyJQ9MnfYtq3lMHUNOVbJFLMiWPYVVbKv7cn9lKWyS4TJlD2yGSGyuT0rqmRzOxXJ01EljTGlY20anpwogLm0jS7YV53z2N/a2gm7LCY9D0OcY6jPtMvC3EPRNcJn23OsMlkdN0nfAv4WeAh4MCI2S3oycCFwIPAt4HURcfcwzTTGrFGK+PSBtcmY8WBt2hHrkzHjoRR9arMQ2T+KiJ+KiM31/qnAFRFxMHBFvW+MMYvG2mSMGSvWJ2NMb3SxSh4NvLzePh+4Cnhnx/Z0tkqmFqBO2SbbWhebw7ZNK12q3bPqbWuPzLFKzhM9aeiFtlP2wVSeruVzrJKpqJLzWCVTdXkB7qUxiDYZY2ZjbcpiNPo0hG2wi/Vu1gLcXdrap2Wzr3MMFUlyEdfaF4tuUyn6lDviFsDnJF0j6eQ6bZ+IuA2g/r33EA00xuzI2iTbNj8FY20yZiRYmx6D9cmYkVCKNuWOuL00Im6VtDdwuaS/yj1BLVYnr5vRGLMuYxeUJdCLNh1wwAFDtc+YDYG1aSrWJ2NGQEn6lNVxi4hb69+3S/o0cBjwPUn7RsRtkvYFbk+U3QJsAZAUjXSmbecswJ0bFbFpOWwb2bBt1MZUeybb27Q7Ns/XTM+JJJlKnycaT+o+d7FHtl00e1ZUyZRVssti3CnbZGoB7qZtcjJfql5HlRyevrRp8+bNvqnGdMTatCNj16chIhCO0ZKXos/rSX3321o8h7pnXert8zNd5vejFH1at0ci6fGSfnxtGzgS+EvgYuCEOtsJwGeGaqQx5lFsR6qwNhkzLqxNj2J9MmZclKJNOSNu+wCfrnvGOwP/JSL+RNLVwKcknQR8B3jtcM00xqwxZkFZMNYmY0aEtWkHrE/GjIhS9GndjltE3Ay8YEr6ncDhfTSirW1ylq0ux+qXEzkxFYWyObTbtMul2j1Zb471sUueeaJN5lgl29oS+1yAe5FWyZwIk7PKpNphq2T/LEKbjDH5WJseZZn6lGtH6zOa4SLLDmHfTOXJbVNfUSUXTU47uizCPquuRd+DUvSpy3IAxpgFM/YhfGPMxsTaZIwZKyXpkztuxqwYpYiPMaYsrE3GmLFSij4treOWske2tU1C2p6WY5tMfZCpsqnh3KZ1btaQb47FcWg75SQpG2rKlpiz3XZx7Fn20pzPtMu5U/bIlJ1yVr6UPXLsVklJRwFnAzsB50TEWRPHdwU+DhwK3AkcGxHfknQEcBawCbgfeHtEXFmXORZ4V13nZyPiHXX6AVQLz+5eHzs1Ii5tnOsA4HrgvRHx271frDFmEEp5MFp1cu1oY7DxzXOutm1axKLUba2WY4m4OUS0yaE+x66Uok/t4twbY5ZO35HbJO0EfBT4eeAQ4HhJh0xkOwm4OyIOAj4MvL9O3w68KiJ+kipC2gV1nXsCHwQOj4jnAvtIWpvX8W7gUxHxQuA44D9MnOvDwB+3vS/GmOXiqJLGmLFSija542bMijHAw9FhwE0RcXNE3A98Ejh6Is/RVKNkABcBh0tSRFwb9VpFwHXAbvXo3LOAGyLijvrY54Fj1i4BeGK9/SRgrTySXg3cXNdljFkh3HEzxoyVUrRpFHPcUpbInEWfIb2odVvbZF9tnUVzaDgnAmSO9bEZATO13Sw7yy7RvI7UwtQ5i0+3jfjY1SrZZSHwtgtz517HEAtwDyQo+wG3NPa3AS9K5YmIByXdA+xJNeK2xjHAtRFxn6SbgGdLOrCu79VUdkqA9wKfk/RW4PHAK+GRtY7eCRwBvK2na1tZ+rKRjPEfUJ8WmTFe30Zk7A87G4muFsAc2+QQ+jTr+zOEra5PG2Pqni3zb2IIy+sY7LXzUJI+jaLjZozJZ07x2UvS1sb+lojYUm9PU9nJk8zMI+m5VPbJI+s23i3pzcCFwMPAl6hG4QCOB86LiA9JeglwgaTnAf8W+HBE/GCMwm+MmU0pD0bGmPIoRZ/ccTNmxZhTfLZHxObEsW3A0xv7+9OwL07k2SZpZyqL410AkvYHPg28MSK+2WjnJcAldZ6TgbVhypOAo+o8X5a0G7AX1SjfayR9gCpwycOS7o2Ij8xzwcaYxVLKg5ExpjxK0adRR5XMWTQb0ra1plUwJ6pk6hw5EYn6tE2mrJKp6JEpe+SmTZum5p9FKpLkfffd98h20zbYjK7YvMfNPG3tjbPK5ESVbH4uOZbIttEzJ/dzbJ19CsYA4nM1cLCkZwLfpQoY8vqJPBdTBR/5MvAa4MqICEm7A58FTouILzYLSNo7Im6XtAfwFuB19aHvUC1Ae56k5wC7AXdExM80yr4X+IE7bcasDqU8GK0680T1y7FE5tjkujCP0yJ17r5sfLOus8s96NLu3OvpsnB2WxvpKrhkStEnBycxZoMTEQ8CpwCXAd+givh4naQzJP2zOtvHgD3ruWu/AZxap58CHAScLumr9c/e9bGzJV0PfBE4KyJuqNN/E/hlSV8D/gA4MUpRVGOMMcaYgbBV0pgVY4g+TlTrqF06kfaexva9wGunlDsTODNR5/GJ9OuBl67Tnveu22hjzKjw+xdjzFgpRZ9G0XFra5WctNWlbIZNO1uOzTA11Ns8X7OeVJsmvxw5X5acaJM5C2037ZHNNs2KKtkkZXdsWiJTi0+3jSo5K7pnKl/KBpljx8zJkxN5cla+lD1y5FEljTGmE9amxdDF9jfL6pdDW2tgn7a/nHM0GXv0w5xzz2NL7HLdQ0eenKd9fX1eJenTKDpuxph8ShEfY0xZWJuMMWOlFH1yx82YFaMU8THGlIW1yRgzVkrRp9FFlcyxRzYtg5AePp1cNHla+ZTNMjU8m2OVnLT9pa6vbZtyFuBO2RubeSbvX84C3H1ZJXMXVU9ZItump75POYtmz7JK5pxjhaJKGmNMZ6xNw9PWGti13iHIjU7ZpU192QFzbX9d6+qLtjbDnOiWOQz1mfZ5n0rRJ0eVNGbFWPNqt/kxxpihsTYZY8bKorVJ0pMlXS7pxvr3Hol8H5B0naRvSPp/tE5v1R03Y1aIeR6M/HBkjBmaZWjTUA9GxpiyWNJz06nAFRFxMHAFjy6j9AiS/g+qKNvPB54H/EPgZbMqdcfNmBXDHTdjzBhZgjYN8mBkjCmPJTw3HQ2cX2+fD7x6WrOA3YBNwK7ALsD3ZlU6iuAkzRuUM98od6X4nDluqXbMmhM2LX9qztS0/fVIzXFrbjfb15z7tssuu0zdzl0OIDUPrDmvrXlfc+a4Ne9TztywWWW6pOfMRcudt5iz5MBQHSZ3xIwxY2QJ2nQ08PJ6+3zgKuCdE3maD0Yi48FozDTvcer5IMWi51W1/T7kho5vS049qTyzrqfL8giLIGe+W19tmieEf1+fby5L0Kd9IuK2+ty3Sdp7Spu+LOlPgduo9OkjEfGNWZWOouNmjMnHHTdjzBiZU5v2krS1sb8lIrZklh3kwcgYUx5z6NO62iTp88BTp5R9V84JJB0EPAfYv066XNLPRsQXUmXccTNmxXDHzRgzRubUpu0RsTl1cBkPRsaY8phDn2ZqU13nK1PHJH1P0r71S6V9gdunZPsF4M8j4gd1mT8GXgyMr+OWGn5OWdXmqbetBa5p78u1Fk4rO/nl6GIXSC1F0NxutjVlm0zZPSdJ2T9ztnNC/acsjZPkhNVPfdY534EuFsrJ/Rz7Z1+dLc9ZM8aMkaG0aRkPRmNmqOUAhqqrC0PY+Jrk2CNn1dO2/DLp8r3Juc7Uc+sse2nO59unXXYJn8vFwAnAWfXvz0zJ8x3glyW9j8oR8DLg38+qNOtpXtLuki6S9Fd1VKaX5EZzMsb0i4OTPIq1yZjxsARtWnswgtkPRi+TtLOkXagejBZilbQ+GTMelvDcdBZwhKQbgSPqfSRtlnROneci4JvAXwBfA74WEZfMqjQ3quTZwJ9ExLOBF1CJ3rrRnIwx/eOO2w5Ym4wZCUvQpkEejHrE+mTMSFj0c1NE3BkRh0fEwfXvu+r0rRHxpnr7oYj4lYh4TkQcEhG/sV6961olJT0R+FngxPok9wP3S8qJ5jTrgtbdTkVjnLyhOba3pp2tGZExxx7ZTM8576x8uUPw09LbRptMWStn0WxfygKYE1Exx5Y4677k5Gu73ZflMjdfzuc+D4V3xLIZSpvGQsmfc8nXtpFZ9OcaEXcCh09J3wo88mAE/MpCG8Zinp0mzrdu/lyrWVsL26qSus7m//GhoiJ2KdvnOdpGI03VmUqfZS9te46ulPCdhbwRt2cBdwC/L+laSedIejwT0ZyAx0RzMsb0j0fcHsHaZMyIsDbtgPXJmBFRijbldNx2Bn4a+N2IeCHwQ1oM7Us6WdLWiZCaxhjTld606Y477hiqjcaYjYn1yRjTOzlRJbcB2yLiK/X+RVTikxPNiajWPNgCIGlqFzbHEplrW2ta91JWwVRExhwrYopZvfO2PfecIe2cyJM5+We1tUukxpx6UuedVVdO+bbb85y3rZWzL8b+JmjB9KZNmzdv9k01pgPWpscwiD4tMhLfUMzTviGuKdfel6JLO9paA4eyvLY9R5d257ZpiO9sSfq07ohbRPwNcIukf1AnHQ5cT140J2NMz9iOVGFtMmZcWJsexfpkzLgoRZty13F7K/CfJW0Cbgb+OVWn71OSTqIKt/vaYZpojGkyZkFZAtYmY0aCtekxWJ+MGQml6FNWxy0ivgpMWz38MdGcckndwJRtclaUnxxLYJftrvRllczJ0+c1dLUTtskzz7mXWbbtOfoUjFLEpw+G0CZjzHxYm3Zk6GenlG2tT3LO15d9c5atrq2FsK3tbx76sjguwj44hNW0T1vnIihFn3JH3IwxI6EU8THGlIW1yRgzVkrRJ3fcjFkhxu69NsZsTKxNxpixUpI+jaLjNnQkPmNKwn8jxpgxYm0anqFsf0Ofbyg74KLvR9tzjMUquMi2juWaJylFn0bRcTPG5FOK+BhjysLaZIwZK6XokztuxqwYpYiPMaYsrE3GmLFSij4tuuO2Hfhh/XujsRcb77o34jVDu+t+RtvKhxAfSUcBZwM7AedExFkTx3cFPg4cCtwJHBsR35J0BHAWsAm4H3h7RFxZlzkWeFdd52cj4h11+gHA+cDu9bFTI+LSWXUNzTXXXLNd0rfZmN/ZjXjNsDGve+W0yTyiT3522jhsxGsG61MWC+24RcRTJG2NiGnhcYtmI173RrxmGPa6h5hgK2kn4KPAEcA24GpJF0fE9Y1sJwF3R8RBko4D3g8cSyWyr4qIWyU9D7gM2E/SnsAHgUMj4g5J50s6PCKuAN4NfCoiflfSIcClwIGpunq92AQR8RTYmN/ZjXjNsDGve9W0yVT42WljXfdGvGawPuXyuGU3wBjTjjUBavOzDocBN0XEzRFxP/BJ4OiJPEdTjZIBXAQcLkkRcW1E3FqnXwfsVo/OPQu4ISLuqI99Hjhm7RKAJ9bbTwJura8rVZcxZgUYQJuMMaYXStEmz3EzZsUYQFD2A25p7G8DXpTKExEPSroH2JMdbQ3HANdGxH2SbgKeLenAur5XU1kgAd4LfE7SW4HHA6+c0qZH6pr/sowxi2TMDzvGmI1NKfq0jI7bliWccwxsxOveiNcMA1/3nOKzl6Stjf0tEbHWzmmxeydPMjOPpOdS2SePrNt4t6Q3AxcCDwNfohqFAzgeOC8iPiTpJcAFkp4XEQ9Pq2vBbMTv7Ea8ZtiY1z1GbTJ5bMTvK2zM696I1wzWpywW3nFrPCxuKDbidW/Ea4bhr3tO8dk+wzu+DXh6Y39/avvilDzbJO1MZXG8C0D6/9u7/9C76jqO489XczpQbKZp6gxXrnQWthBLhP5wptNKM5W2Iq1GQqYoCbalRYiCFqUSGdi0RMQplrilpfkDArOlNsvmms4VurK2kVYEKdNXf5zP3E237+62773fcz/n9YDL995zPufXzuW1z/ueX5oG3AGcafuZnvVcCiwtbc4GXimj5gNzSpuHJU2huSh53dbmNSxd/M52cZuhm9vd0myKPnTx+wrd3O4ubjMkn/qVa9wiRsiOXEPSR1g9AsyQNF3SrsBcYMnr2iwBzirvTwcesG1JU4G7gIW2H+qdQNK+5e9ewDnAojLqWWB2GXcYMAVYP9a8IqLdBpRNERE7raZsSuEW0XG2NwLn0tzFcSXNHR9XSLpU0sml2fXA3uXatS8DC8rwc4FDgK9Jery89i3jrpH0JPAQcIXtp8rwC4EvSPodcAvwWTcpOda8IiIiIjptqIWbpDmSVklaLWnBtqcYPZIOkvSgpJWSVkg6vwx/i6RfSHq6/N1rotd1ECRNkrRc0k/L5+mSlpXtvrUc0amGpKmSbpf0x7LPjx70vh7EL0e277b9LtvvtH15GfZ120vK+//aPsP2IbaPsr2mDL/M9u6239fzWlfGzbM9s7wW9yzrSdvH2D6itL93W/MatC5kE3Q7n7qWTTD8fKrpV+026UI+JZuSTW3rO7XV0Ao3bX5W1InATGCemmc41WYjcKHtw4APAl8q27kAuN/2DOB+Nh+xqM35NEdtNrkSuKps9ws01zfV5Brg57YPBY6g2faB7ut0jsZXh7IJup1PXcsmGHI+JZvGX4fyKdm0WbKpBX2nthrmEbd+nhU18mw/b/u35f2/ab6MB/L/z8G6keb26FUpN5b4COVaJkkCjqV57hdUtt2S9gQ+RHMaIbZftv0iA97X6RyNu05kE3Q3n7qWTTAx+ZRsGohO5FOyKdlEy/pObTXMwm1Lz4o6cIjLHzo1z7CaBSwD9rP9PDQBBdR47c7VwEU0t3+H5jlfL7q5hgrq2+fvANYDPyynOSyStDsD3tfpHI27zmUTdC6fupZNMAH5lGwaiM7lU7Ip2dSGvlNbDbNw6+dZUdWQtAfwY+AC2/+a6PUZNEkfBdbZfqx38Baa1rTPdwHeD3zf9izgPwz4NI4d6Ri1OYBaovbv6Rt0KZ86mk0w5HxKNg1MF76rr0k2dWJ/j0Tfqa2GWbj186yoKkiaTBM8N9v+SRn8d0n7l/H7A0O56cIQHQOcLOnPNKdyHEvzS9JUNc/9gvr2+Vpgre1l5fPtNGE00H1dUwC1RGeyCTqZT13MJpiAfEo2DURn8inZlGxqU9+prYZZuPXzrKiRV85Pvh5Yafs7PaN6n4N1FnDnsNdtkGwvtD3N9sE0+/YB258GHqR57hdUtt22/wY8J+ndZdBs4EkGvK/TORp3ncgm6GY+dTGbYGLyKdk0EJ3Ip2RTsomW9Z3aapdtNxkftjdK2vSsqEnADbZXDGv5Q3QM8BngCUmPl2FfBa4AbpM0n+YBxGdM0PoN21eAxZIuA5ZTLkatyHnAzeU/1DXA52h+EBnYvm5zoIyiDmUTJJ961Z5NMOR8SjaNvw7lU7Jps2RT+k5bpVo2JKILJk+e7KlTp273dBs2bHjM9pEDWKWIiGRTRLTWjuRTW7NpaEfcImLntf0QfkR0U7IpItqqpnxK4RYxYmoJn4ioS7IpItqqlnxK4RYxYmoJn4ioS7IpItqqlnxK4RYxYmoJn4ioS7IpItqqlnxK4RYxYmoJn4ioS7IpItqqlnxK4RYxQmq6wDYi6pFsioi2qimfUrhFjJhawici6pJsioi2qiWfUrhFjJhawici6pJsioi2qiWf3jTRKxARERERERFjyxG3iBFTy69GEVGXZFNEtFUt+ZTCLWLE1BI+EVGXZFNEtFUt+ZTCLWKE1HRnpIioR7IpItqqpnxK4RYxYmoJn4ioS7IpItqqlnxK4RYxYmoJn4ioS7IpItqqlnxK4RYxYmoJn4ioS7IpItqqlnxK4RYxYmoJn4ioS7IpItqqlnzKc9wiRsimC2y39xURMUjJpohoq4nIJklnSFoh6VVJR47Rbo6kVZJWS1qwrfmmcIsYMYMIoG0Fh6TdJN1axi+TdHAZ/mFJj0l6ovw9tmeaT0r6fQmub/YMf7ukByUtL+NP6hm3sCxjlaQTdvKfKiKGaNido0F1jCKiPhPwo9IfgE8Av9xaA0mTgO8BJwIzgXmSZo410xRuESNmvDtHfQbHfOAF24cAVwFXluEbgI/Zfi9wFnBTmefewLeA2bYPB/aTNLtMcwlwm+1ZwFzg2jLNzPL5cGAOcG1Zt4gYARPwq/ZAOkYRUZ9hF262V9petY1mRwGrba+x/TKwGDhlrAlSuEWMmAF0jvoJjlOAG8v724HZkmR7ue2/luErgCmSdgPeATxle30Zdx9w2qZNAPYs798MbJr+FGCx7Zds/wlYXdYtIkbAsAu3QXWMIqI+E3DErR8HAs/1fF5bhm1Vbk4SMWIGEChbCo4PbK2N7Y2S/gnsTXPEbZPTgOW2X5K0Gji0nFK5Fvg4sGtp9w3gXknnAbsDx/Us49evW48xAywi2mOInZ3t0U++RUTldiCf9pH0aM/n62xf19tA0n3A27Yw7cW27+xjGdrCsDFXNIVbxGi5B9hnB6abMkYA9RMcY7aRdDjN6ZPHA9h+QdIXgVuBV4Ff0RyFA5gH/Mj2tyUdDdwk6T19rkdEtNMgsmlCOkYRUZ0dyacNtueM1cD2cWON78Na4KCez9PYfBbSFqVwixgh2wqRHdRPcGxqs1bSLjSnOP4DQNI04A7gTNvP9KzrUmBpaXM28EoZNZ/mGjZsPyxpCk2gbneARUQ7DCibJqRjFBF1GVQ+jYNHgBmSpgN/obnO/1NjTZBr3CLiteCQtCtNcCx5XZslNDcfATgdeMC2JU0F7gIW2n6odwJJ+5a/ewHnAIvKqGeB2WXcYcAUYH1ZxtxyB8vpwAzgN+O6pRHRNf3kW0TEuJJ0qqS1wNHAXZLuKcMPkHQ3NJeeAOfSHBFcSXPjthVjzrel56RHxBCVW/JfDUwCbrB9uaRLgUdtLylHxW4CZtEcaZtre42kS4CFwNM9szve9jpJtwBHlGGX2l5cljUT+AGwB80pSxfZvreMuxj4PLARuMD2zwa75RExqiSdCnwXeCvwIvC47RMkHQAssn1SafeGfJuodY6I2Bkp3CIiIiIiIloup0pGRERERES0XAq3iIiIiIiIlkvhFhERERER0XIp3CIiIiIiIlouhVtERERERETLpXCLiIiIiIhouRRuERERERERLZfCLSIiIiIiouX+B9sz5+SjPqj0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x360 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(15,5))\n",
    "plt.subplot(1,3,1)\n",
    "im = plt.imshow(U, cmap='gray')\n",
    "plt.colorbar(im)\n",
    "plt.title('finite diffenence u')\n",
    "\n",
    "plt.subplot(1,3,2)\n",
    "im = plt.imshow(F, cmap='gray')\n",
    "plt.colorbar(im)\n",
    "plt.title('true f')\n",
    "\n",
    "plt.subplot(1,3,3)\n",
    "im = plt.imshow(conv_U, cmap='gray')\n",
    "plt.colorbar(im)\n",
    "plt.title('conv u')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 246,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.1200720455343632e-06"
      ]
     },
     "execution_count": 246,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean((conv_U - F)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "0.00011501938590870051"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 217,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([   0.     , -120.85266,    0.     ,    0.     , -362.55798,\n",
       "          0.     , -241.70532,    0.     ,    0.     ,  120.85266,\n",
       "       -120.85266,    0.     , -120.85266,  241.70532,    0.     ,\n",
       "       -241.70532, -120.85266,    0.     ,  241.70532, -241.70532,\n",
       "        120.85266,    0.     ,  120.85266, -120.85266,  120.85266,\n",
       "        120.85266,  120.85266,    0.     ,    0.     ,  120.85266,\n",
       "       -120.85266,    0.     ,  120.85266,    0.     , -120.85266,\n",
       "        120.85266, -241.70532,  120.85266,    0.     ,    0.     ,\n",
       "       -120.85266, -241.70532,  362.55798, -120.85266, -120.85266,\n",
       "        120.85266, -241.70532, -241.70532,  120.85266, -241.70532,\n",
       "        241.70532, -120.85266,  241.70532, -362.55798,  241.70532,\n",
       "        241.70532, -241.70532,  241.70532, -241.70532,    0.     ,\n",
       "          0.     , -120.85266,  120.85266,    0.     ], dtype=float32)"
      ]
     },
     "execution_count": 217,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "conv_U[1,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 509,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.51223058917094e-08"
      ]
     },
     "execution_count": 509,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(U[1:4, 3:6] * np.array([[0,  1,  0],\n",
    "                        [1, -4,  1],\n",
    "                        [0,  1,  0]]))/h**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 394,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = torch.arange(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 395,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[0, 1, 2]])"
      ]
     },
     "execution_count": 395,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.unsqueeze(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 304,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 2])"
      ]
     },
     "execution_count": 304,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 305,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  0, 111,   2])"
      ]
     },
     "execution_count": 305,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
